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ABSTRACT 


HARMS,  DEWEY  ELVIN.  Four-Dimensional  Data  Assimilation  of  GALE  Data  using  a 
Multivariate  Analysis  Scheme  and  a  Mesoscale  Model  with  Diabatic  Initialization.  (Under 
the  direction  of  Sethu  Raman). 

A  method  of  assimilating  3-hourly  sounding  data  is  developed  and  successfully 
tested  in  this  study.  First,  the  successive  corrections  scheme  of  Bratseth  (1986),  which 
converges  to  optimum  interpolation,  is  applied  for  the  numerical  analysis  of  data  collected 
during  the  Genesis  of  Atlantic  Lows  Experiment  (GALE).  Univariate  analyses  of  the 
mass  and  wind  field  are  produced.  The  coupling  of  the  mass  and  wind  field  is  achieved 
by  further  iterations  of  the  geopotential  utilizing  improving  estimates  of  the  geostrophic 
wind  to  extrapolate  the  geopotential  to  the  grid  points.  The  univariate  wind  analysis  is 
then  corrected  for  the  new  geostrophic  wind. 

Next,  diabatic  forcing  is  incorporated  into  a  vertical  mode  initialization  scheme  to 
provide  more  realistic  initial  conditions  and  to  shorten  the  spinup  time  of  the  Naval 
Research  Laboratory  /  North  Carolina  State  University  (NRL/NCSU)  mesoscale  model. 
Latent-heating  profiles  are  computed  from  'spun-up'  model-generated  and  observed 
rainfall.  The  latent  heating  is  distributed  in  the  vertical  according  to  the  cumulus 
convective  parameterization  scheme  (Kuo  scheme)  of  the  model.  Compatibility  between 
the  specified  heating  during  initialization  and  the  heating  during  early  model  integration  is 
retained  by  merging  the  model  integrated  rainfall  and  heating  rates  with  those  rates  from 
the  initialization. 

Finally,  the  multivariate,  successive  correction  analysis  scheme  and  the  diabatic 
initialization  procedure  are  combined  with  the  NRL/NCSU  model  to  form  an  intermittent 
data-assimilation  system.  Assimilations  of  the  GALE  data  over  a  2  1/2-day  period  were 
performed  with  differing  update  cycles  of  3,  6,  and  12  h.  Twelve-hour  NMC 


hemispheric  analyses  served  as  the  "no  assimilation"  control  case  for  comparison.  The 
assimilation  of  3-hourly  GALE  data  led  to  large  decreases  in  background  forecast  rms 
errors  and  smaller  decreases  in  analysis  rms  error.  Better  consistency  in  time  was 
achieved  between  forecasts  and  analyses  in  the  assimilation  experiments.  Rainfall 
prognoses  from  the  assimilated  states  verified  reasonably  well  with  the  observed  rainfall 
and  showed  much  more  rapid  spinup  and  better  overall  patterns  than  did  the  "no 
assimilation"  precipitation  forecasts. 
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1.  INTRODUCTION 


Four-dimensional  data  assimilation,  in  the  context  of  numerical  weather 
prediction,  is  the  process  of  specifying  initial  conditions  for  the  integration  of  a  numerical 
forecast  model.  In  this  process,  new  observational  data  are  merged,  using  some  analysis 
procedure,  with  the  ongoing  integration  of  the  prediction  model.  Data  analysis  and  model 
prediction,  along  with  model  initialization,  make  up  the  major  components  of  the  data- 
assimilation  process.  A  short-range  model  forecast  serves  as  a  first  guess  for  the 
subsequent  analysis. 

The  overall  purpose  of  this  research  is  to  develop  a  method  of  data  assimilation, 
which  includes  a  cost-effective  multivariate  analysis  scheme,  accounts  for  diabatic  effects 
in  the  model  initialization,  and  uses  a  fine-resolution  mesoscale  model,  in  order  to 
produce  improved  analyses  of  mesoscale  circulations  and  then  more  reliable  forecasts  of 
coastal  frontogenesis/cyclogenesis  and  associated  precipitation.  The  goal  here  is  also  to 
investigate  the  impact  of  assimilating  high-frequency  observational  data  using  a 
mesoscale  model.  To  achieve  this  goal,  the  following  objectives  were  required:  1) 
develop  and  test  a  computationally  efficient,  multivariate  objective  analysis  scheme,  2) 
shorten  model  spinup  by  including  diabatic  heating  in  the  initialization  procedure,  and  3) 
combine  these  components  with  a  mesoscale  numerical  model  to  form  an  intermittent 
data-assimilation  system  and  then  evaluate  the  results. 

Four-dimensional  data-assimilation  methods,  along  with  the  most  commonly  used 
objective  analysis  and  initialization  techniques,  are  examined  from  a  historical  perspective 
in  Chapter  2.  Operational  techniques,  including  intermittent  data  assimilation  and 
Newtonian  nudging,  and  next-generation  methods  (Kalman-Bucy  filtering  and  the  adjoint 
method)  are  briefly  described.  Several  methods  are  compared,  with  primary  emphasis 


2 


being  placed  on  recent  papers  dealing  with  the  operational  assimilation  techniques. 
Ongoing  and  future  research  is  outlined,  and  some  important  implications  of  this  research 
are  discussed.  Following  this  extensive  literature  review,  the  NRL/NCSU  Limited-Area 
Numerical  Weather  Prediction  Model  (refeired  to  hereafter  as  the  NRL/NCSU  mesoscale 
model)  is  described  in  Chapter  3.  This  model,  which  was  developed  jointly  by  the  U.S. 
Naval  Research  Laboratory  (NRL)  and  North  Carolina  State  University  (NCSU),  is  used 
throughout  the  course  of  this  research.  The  description  includes  the  model's  dynamics, 
numerics,  physics,  and  vertical  mode  initialization  procedure. 

In  Chapter  4,  the  successive  correction  scheme  of  Bratseth  (1986)  is  applied  for 
the  multivariate  analysis  of  data  colleaed  during  the  Genesis  of  Atlantic  Lows  Experiment 
(GALE),  which  was  conducted  over  the  southeastern  U.S.  during  the  winter  of  1986. 
The  newly-developed  multivariate  analysis  interpolation  procedure  is  described  in  detail, 
and  the  utility  of  the  analysis  scheme  is  measured.  After  univariate  analyses  of  the  mass 
and  wind  fields  are  produced,  the  coupling  of  these  fields  is  achieved  through  additional 
iterations  of  the  geopotential,  utilizing  improved  estimates  of  the  geostrophic  wind  to 
extrapolate  the  geopotential  to  the  grid  points.  The  wind  analysis  is  then  corrected  for  the 
new  geostrophic  wind.  This  study  revealed  that  the  objective  analysis  scheme 
successfully  resolved  mesoscale  structures  and  signiHcantly  improved  the  absolute 
vorticity  and  vertical  velocity  fields. 

Diabatic  forcing  from  convective  heating,  as  applied  in  the  model  initialization 
procedure,  is  discussed  in  Chapter  5.  Latent-heating  profiles  are  computed  from  'spun- 
up'  model-generated  and  observed  rainfall.  The  latent  heating  is  distributed  in  the  vertical 
according  to  the  cumulus  convective  parameterization  scheme  (Kuo  scheme)  of  the 
model.  Compatibility  between  the  specified  heating  during  initialization  and  the  heating 
during  early  model  integration  is  retained  by  merging  the  model  integrated  rainfall  and 
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heating  rates  with  those  rates  from  the  initialization.  The  inclusion  of  convective  heating 
in  the  model  initialization  substantially  reduced  the  spinup  of  model  rainfall. 

In  Chapter  6,  the  multivariate,  successive  correction  analysis  and  the  diabatic 
initialization  components  are  combined  with  a  mesoscale  numerical  weather  prediction 
model  to  form  an  intermittent  data-assimilation  system.  A  thorough  evaluation  of  the 
assimilation  system  was  performed  using  GALE  data.  Finally,  a  summary  and  major 
conclusions  are  presented  in  Chapter  7. 
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2.  LITERATURE  REVIEW  OF  DATA  ASSIMILATION 

2.1  Introduction 

Two  major  motivations  for  using  data  assimilation  exist:  as  an  analysis/diagnostic/ 
research  tool  and  for  operational  weather  forecasting.  Data  assimilation  has  been  applied 
not  only  in  meteorology  (air  pollution  and  planetary  boundary  layer  studies,  forecast  case 
studies,  quantitative  assessment  of  new  observing  systems,  among  others),  but  also  in 
oceanography  for  describing  ocean  currents  (Ghil  1989;  Robinson  1986).  To  summarize 
all  relevant  work  pertaining  to  data  assimilation  is  indeed  a  difficult  task  because  the 
contents  are  spread  so  widely.  The  purpose  of  this  review  is  only  to  provide  an  overview 
of  four-dimensional  data  assimilation  with  primary  emphasis  on  assimilation  methods 
currently  useful  for  operational  weather  fraecasting.  To  make  this  review  more  complete, 
attention  has  also  been  given  to  the  state-of-the-art  or  "next-generation"  techniques,  but  to 
a  lesser  extent.  Excellent  reviews  covering  the  early  days  of  data  assimilation  research  in 
both  simulation  and  real-data  studies  are  available  in  the  literature.  The  interested  reader 
is  referred  to  Bengtsson  (1975a),  McPherson  (1975),  Hollingsworth  (1986),  and  Bourke 
et  al.  (1985).  Plus,  a  comprehensive  history  of  data  analysis  and  assimilation  is  given  by 
Daley  (1991),  whereas  a  more  abbreviated  examination  of  data  assimilation  techniques  is 
rendered  by  Harms  et  al.  (1992b). 

Numerical  weather  prediction  (NWP)  has  classically  been  viewed  as  an  initial- 
value  problem  where  the  governing  equations  of  geophysical  fluid  dynamics  are 
integrated  forward  in  time  from  a  set  of  initial  values.  The  quality  of  NWP  is  strongly 
dependent  on  the  accuracy  of  specifying  these  initial  conditions  and  on  the  ability  to 
model  mathematically  the  dynamics  and  physical  processes  of  the  atmosphere.  In  a 
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pioneering  paper,  Chamey  et  al.  (1969)  suggested  combining  past  and  current  data  in  a 
numerical  model  such  that  the  model's  equations  provide  time  continuity  and  dynamic 
coupling  among  the  atmospheric  fields.  This  concept,  which  has  merged  objective 
numerical  analysis  and  numerical  weather  prediction,  has  become  known  as  four¬ 
dimensional  data  assimilation  (FDDA)  and  has  proven  to  be  a  major  advance  in  NWP 
during  the  past  20  years. 

What  inspired  this  major  advance,  or  in  other  words,  what  spurred  this 
development  of  data  assimilation?  The  advent  of  meteorological  satellites  in  the  1960s 
raised  the  possibility  that  nearly  continuous  atmospheric  temperature  observations  would 
be  available  on  a  global  basis.  However,  these  space-based  observing  systems  measure 
only  radiance  (or  temperature)  distributed  in  space  and  time,  rather  than  at  fixed  locations 
and  times.  In  order  to  fully  utilize  this  new  source  of  data,  the  numerical  weather 
prediction  techniques  had  to  be  adjusted.  Beginning  with  Chamey  et  al.  (1969), 
Smagorinsky  et  al.  (1970),  Rutherford  (1972),  and  Morel  and  Talagrand  (1974), 
research  progressed  in  reconstructing  unobserved  variables  from  the  observed  variables 
through  the  numerical  model's  dynamical  coupling  between  those  variables.  By 
combining  information  about  the  state  of  the  atmosphere,  earlier  observations  are  carried 
forward  to  provide  an  independent  source  of  information  to  be  added  to  the  newly 
acquired  observations. 

Morel  (1981)  further  illustrated  why  data  assimilation  is  essential  in  NWP  by 
listing  five  key  reasons:  1)  the  inadequate  distribution  (spatial  gaps)  of  the  twice-daily 
upper-air  sounding  data,  2)  the  discrepancy  between  the  conventional  observations  as 
point  measurements  and  the  true  volume  averages  required  by  numerical  models,  3)  the 
inherently  asynoptic  character  of  remote  observations  obtained  from  sunsynchronous 
orbiting  satellites,  4)  the  inadequate  vertical  resolution  of  remote  observations  of  cloud 
motions  from  geostationary  satellites,  and  5)  the  significant  random  and  systematic 
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em)rs  involved  in  the  data  processing  required  for  reconstructing  atmospheric  fields  from 
remotely  measured  physical  quantities.  He  concluded  that  any  weather  prediction  model 
must  be  initialized  by  merging  the  new  observations  with  the  currently  estimated 
meteorological  fields,  computed  on  the  basis  of  earlier  observations,  while  taking  into 
account  the  dynamical  constraints  between  successive  model  states,  specified  by  the 
governing  dynamical  equations. 

This  process  of  merging  new  observational  data  with  the  ongoing  integration  of  a 
numerical  forecasting  model  is  known  as  "data  assimilation"  or,  equivalently,  "four- 
dimensicnal  data  assimilation"  in  consideration  of  the  time-space  distribution  of  the  data 
base.  Then,  as  stated  by  Warner  (1987),  "the  overall  objective  of  the  assimilation 
process  is  to  provide  the  best  possible  initial  state  from  which  to  begin  a  forecast,  where 
the  term  "best"  implies  an  appropriate  balance  as  well  as  reflects  optimal  use  of  four¬ 
dimensional  data  (synoptic  as  well  as  asynoptic)  to  define  the  structures  on  all  scales  at 
the  initial  time." 

In  the  past  20  years,  many  techniques  have  been  developed  to  insen  data  into 
dynamic  models.  However,  the  first  data-assimilation  experiments  were  simplistic  and  of 
limited  success.  The  observation  simply  replaced  the  model  forecast  value  at  the  model 
grid  point  nearest  the  observation  location  (Jastrow  and  Halem  1970).  This  technique, 
known  as  direct  insertion,  is  inadequate  from  the  initialization  point  of  view  (Bengtsson 
1975b).  Direct  insertion  shocks  the  system  and  generates  large-amplitude  waves,  or 
gravitational  oscillations.  For  example,  if  observations  of  the  mass  field  are  inserted  into 
a  primitive  equation  model,  an  imbalance  is  created  between  the  mass  and  wind  fields. 
When  the  model  integration  is  resumed,  this  imbalance  is  manifested  as  large-amplitude 
gravitational  oscillations;  this  is  the  model's  attempt  to  restore  the  dynamic  balance  which 
was  disrupted  by  the  data  insertion.  Techniques  such  as  damping  time  integration 
schemes  and  time  filters  were  developed  to  dampen  these  non-meteorological  waves. 
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However,  this  damping  must  be  rather  strong  and  can  be  harmful  to  the  meteorological 
modes.  Bengtsson  (1975b)  stated  that  the  shock  effect  will  be  reduced  if  locally  analyzed 
data  (interpolating  the  observation  to  several  nearby  grid  points)  are  insened  into  the 
model  instead  of  observations  only  (indirect  insertion).  For  example,  applying  a  local 
multivariate  analysis,  or  correspondingly,  relating  wind  and  geopotential  by  the 
geostrophic  relation  and  inserting  both  height  and  winds  simultaneously  proved 
successful  by  speeding  the  updating  considerably  (see  Fig.  2.1).  The  multivariate 
procedure  mentioned  here  produces  a  simultaneous  weighting  of  mass  and  motion 
observations,  subject  to  the  constraint  of  geostrophy  (Petersen  1968;  Eddy  1973;  Daley 
and  Puri  1980). 

Of  the  several  FDDA  methods  that  have  been  investigated  over  the  last  two 
decades,  some  have  been  implemented  operationally  while  others  have  not  yet  been  used 
or  are  in  the  developmental  stage.  In  the  1970s,  regionaVmesoscale  weather  forecasting 
models  were  developed,  normal  mode  initialization  was  introduced,  and  global  data 
assimilation  became  operational.  With  the  advent  of  mesoscale  modeling,  FDDA 
applications  are  also  geared  toward  the  mesoscale,  with  the  emphasis  shifting  from 
applications  on  the  global  scale. 

In  operational  numerical  prediction  systems,  FDDA  can  be  categorized  into  two 
broad  areas:  continuous  methods  and  intermittent  (analysis-forecast  cycle)  methods.  The 
former  refers  to  the  insertion  of  data  into  a  model  as  it  is  received,  in  a  temporally 
continuous  fashion.  The  analysis-forecast  cycle  clearly  illustrates  the  four  components  of 
data  assimilation:  quality  control,  objective  analysis,  initialization,  and  an  initial  guess 
from  a  short-range  forecast.  With  this  method,  the  data  are  assimilated  intermittently  at 
specified  intervals.  So-called  next-generation  data-assimilation  methods  are  being 
researched  today,  with  the  two  most  prevalent  being  the  adjoint  method  based  on  the 
variations  of  calculus  and  the  Kalman-Bucy  filter  technique. 
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Figure  2.1.  Successive  reduction  of  wind  error  using  a  direct  insertion  (upper  curve)  of 
hei^t  data  in  a  batotropic  model.  Conesponding  results  using  height  as  weU  as  from  the 
derived  geostrophic  winds  Oower  curve)  (McPherson  1975). 
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Since  quality  control,  objective  analysis,  and  initialization  are  intricate  parts  of  data 
assimilation,  it  is  appropriate  to  consider  each  of  these  components  in  the  present 
discussion.  Quality-control  methods  and  the  most  commonly  used  objective  analysis  and 
initialization  techniques  are  overviewed,  followed  by  a  discussion  of  various  data- 
assimilation  methods.  Then,  two  studies,  comparing  operational  data-assimilation 
techniques,  are  reviewed,  and  some  important  implications  of  FDDA  and  ongoing  and 
future  research  in  data  assimilation  are  outlined. 

2.2  Quality  Control 

Quality  control  is  an  integral  part  of  a  data-assimilation  system.  Quality-control 
algorithms  are  designed  to  modify  or  reject  erroneous  meteorological  data.  Following 
Daley  (1991),  observational  data  errors  can  be  classified  into  two  groups:  natural  error 
and  gross  (or  rough)  error.  The  natural  error  includes  instrument  error  and  error  of 
representativeness.  The  data  describe  the  behavior  of  the  instrument  itself,  not  the 
behavior  of  the  meteorological  parameter  it  is  intended  to  measure.  Every  instrument  is 
approximate  by  its  very  nature.  Errors  of  representativeness  are  deviations  caused  by 
small-scale  perturbations  and  are  also  referred  to  as  micrometeorological  errors.  Gross 
errors  originate  from  improperly  calibrated  or  malfunctioning  instruments,  incorrect 
registration  of  observations,  incorrect  coding  of  observations,  and  telecommunication 
errors.  These  errors  (natural  and  gross)  can  be  either  random  or  spatially  or  temporally 
correlated,  and  there  can  be  systematic  biases.  The  systematic  errors  can  result  from 
improper  calibration  of  an  instrument  or  from  the  influence  of  some  persistent  factor  that 
is  not  accounted  for  or  is  accounted  for  inaccurately. 

Several  quality-control  techniques  are  used  routinely  to  screen  bad  data.  These 
techniques  can  be  divided  into  four  major  categories  (following  Gandin  1988).  The  first 
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two,  plausibility  check  and  check  for  contradictions,  are  used  to  identify  gross  errors 
based  on  the  physical  reasonableness  of  the  data.  Plausibility  checks  are  the  most  widely 
used  quality-control  methods.  These  checks  analyze  each  datum  independently  of  other 
data.  A  simple  plausibility  check  is  one  that  rejects  data  values  that  can  never  occur  in 
reality— for  example,  positive  temperatures  (*C)  at  300  mb.  Other  versions  compare  the 
datum  with  the  climatological  mean  or  with  a  background  field  (numerical  forecast).  If 
the  deviation  is  too  large,  the  datum  may  be  rejected.  The  check  for  contradictions  is 
based  on  an  analysis  of  two  or  more  parameters  at  the  same  point.  An  example  is  the 
occurrence  of  rain  in  the  absence  of  clouds. 

The  remaining  quality-control  procedures  include  checks  that  rely  on  some 
common  information  and  redundancy  between  observations.  The  spatial  continuity  (or 
consistency)  check  compares  a  datum  with  data  at  adjacent  locations.  This  method,  also 
known  as  the  buddy  check,  demands  spatial  consistency  among  the  data.  A  temporal 
continuity  check  can  also  be  made  in  which  temporal  continuity  is  required  with  past 
observations.  Finally,  checks  using  dynamic  relations  such  as  the  geostrophic  or 
hydrostatic  relation  can  be  used  to  check  geopotential  height  with  temperature  and  mass 
with  wind.  These  checks  require  the  data  to  obey  the  dynamic  relation  at  least 
approximately,  otherwise  the  data  are  suspected  or  rejected. 

The  necessity  for  including  automated  quality  control  in  operational  NWP  was 
recognized  during  the  early  era  of  NWP.  Significant  advances  have  been  made  in  quality 
control,  including  the  development  of  methods  using  optimum  interpolation  (Rutherford 
1976;  Lorenc  et  al.  1977;  Lorenc  1981).  Recent  advances  include  the  complex  quality- 
control  method  of  Gandin  (1988)  and  the  Bayesian  approach  of  Lorenc  and  Hammon 
(1988).  Reviews  on  various  quality  control  methods  have  been  authored  by  Belousov  et 
al.  (1968),  Gustavsson  (1981),  and  Lorenc  (1985). 


2.3  Objective  Analysis  Techniques 


Objective  analysis  (the  second  major  component  of  FDDA)  and  quality  control 
have  become  intertwined  and,  as  a  result,  the  data-assimilation  process  is  more  internally 
consistent.  Excellent  reviews  on  objective  analysis  have  been  presented  by  McPherson 
(1976)  and  Gustavsson  (1981).  The  present  discussion  contains  a  synopsis  of  selected 
topics  from  these  two  reviews.  According  to  McPherson  (1976),  objective  analysis  is  a 
process  in  which  meteorological  observations  distributed  in  space  and  time  are  combined 
with  forecasts  from  previous  analyses  and  perhaps  with  climatology  to  form  a  numerical 
representation  of  the  state  of  the  atmosphere.  This  representation  takes  the  form  of  values 
of  the  meteorological  variables  at  regularly  spaced  grid  points  subject  to  various 
mathematical  and  physical  constraints. 

The  objective  analysis  process  consists  of  three  subprocesses  which  are  essential 
for  the  overall  success  or  effectiveness  of  the  FDDA  system:  1)  filtering  of  small 
amplitude  random  and  systematic  errors;  2)  interpolation  to  a  regular  network  of  grid 
points  or,  in  the  case  of  spectral  analysis,  integration  of  the  representing  mathematical 
space  functions  over  the  irregularly  spaced  observations;  and  3)  forced  adjustment  of  the 
meteorological  variables  using  dynamic  relationships  among  these  variables. 

Within  a  relatively  short  period  of  time,  objective  (numerical)  analysis  schemes 
were  independently  developed  by  several  meteorologists.  Panofsky  (1949)  devised  the 
first  objective  analysis  method-polynomial  inteipolation,  or  the  so-called  surface-fitting 
type.  An  extension  of  this  method  was  developed  by  Gilchrist  and  Cressman  (1954)  and 
became  the  flrst  operational  objective  analysis.  In  this  technique,  mathematical 
(polynomial)  functions  are  adjusted  to  observed  data  in  the  close  vicinity  of  the  grid 
point,  for  which  analyzed  values  are  required.  The  adjustment  or  fit  is  obtained  by  a  least 
squares  technique.  The  polynomial  method  is  nonlinear  since  nonlinear  functions  are 
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used  to  approximate  the  variation  of  the  analyzed  variable.  However,  the  resulting 
analyzed  value  at  the  grid  point  is  a  linear  function  of  the  observed  data  and,  in  this 
respect,  this  method  is  similar  to  the  other  analysis  methods  that  will  be  described. 

Bergthorsson  and  Doos  (1955)  introduced  the  successive  correction  method,  and  a 
similar  method  was  devised  by  Cressman  (1959).  Qessman's  successive  correction 
technique  essentially  replaced  the  polynomial  interpolation  method  because  the  latter 
produced  unreasonable  analyses  at  the  edge  of  data-rich  and  data-sparse  areas.  In  the 
successive  correction  technique,  a  forecast  model  provides  the  preliminary  estimate  (first 
guess)  of  the  field  to  be  analyzed.  The  basic  idea  of  the  method  is  to  correct  this 
preliminary  field  iteratively  during  several  analysis  "scans";  the  results  of  one  scan 
become  the  first  guess  for  the  next  scan.  The  estimate  is  modified  by  a  combination  of 
corrections  computed  for  each  grid  point  The  corrections,  which  are  proportional  to  the 
difference  between  the  observed  and  first  guess  values,  are  weighted  empirically  with 
observations  nearest  the  grid  point  weighted  the  most  heavily. 

By  using  a  short-range  forecast  based  on  a  previous  analysis  as  the  first  guess,  the 
effect  of  past  data  is  allowed  to  influence  the  analysis  and  thereby  contributes  to  the 
temporal  continuity  and  spatial  coherence  of  the  analysis.  For  many  years,  successive 
correction  methods  were  most  widely  used  in  operational  forecast  centers. 

This  early  dominance  of  the  successive  correction  methods  was  taken  over  during 
the  last  decade  by  statistical  (optimal)  interpolation  schemes,  which  were  originally 
introduced  by  Eliassen  (1954)  and  Gandin  (1963).  Statistical  interpolation  was  also 
studied  extensively  by  Alaka  and  Elvander  (1972)  and  Phillips  (1976).  This  analysis 
technique  is  based  on  statistical  linear  regression  and  provides  a  systematic  framework 
for  blending  observations  of  differing  error  characteristics  with  recent  predictions  or 
climatology.  More  accurate  data  receive  more  weight  in  the  analysis.  As  in  the 
successive  correction  method,  the  analyzed  value  at  the  grid  point  is  the  sum  of  the  first 
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guess  and  a  linear  combination  of  corrections,  which  are  proportional  to  the  difference 
between  observational  and  first  guess  values.  The  weighting  coefficients  are  determined 
from  the  condition  that  the  mean- square-error  of  the  analyzed  values  be  minimum,  and 
they  depend  on  the  spatial  covariances  among  the  analyzed  variables.  This  method  is,  in 
principle,  spatially  coherent  and,  like  its  counterpart  of  successive  correction, 
incorporates  temporal  continuity  through  the  use  of  a  short-range  forecast  from  the 
preceding  analysis  as  a  first  guess. 

The  polynomial  interpolation  method  mentioned  earlier  reappeared  in  the  British 
Meteorological  Office  analysis  (Dixon  et  al.  1972)  and  in  the  spectral  analysis  method 
devised  by  Flattery  (1971).  In  the  spectral  analysis  technique,  mathematical  functions  are 
globally  adjusted  to  fit  observed  data.  One  additional  analysis  method  that  has  been 
utilized  in  the  FDDA  context  (see  discussion  in  section  5)  is  the  variational  technique 
introduced  by  Sasaki  (1958).  This  method  is  a  post-analysis  adjustment  technique  based 
on  the  calculus  of  variations  and  is  very  effective  in  making  the  analyzed  fields 
compatible  with  a  forecast  model.  A  functional  is  used,  which  minimizes  the  analyzed- 
minus-observed  difference,  filters  undesirable  high-frequency  and  high-wavenumber 
features,  and  employs  dynanucal  constraints.  These  constraints  may  be  strong  (satisfied 
exactly)  or  weak  (satisfied  {q)proximately). 

Over  the  past  30  years,  numerous  variations  of  these  different  objective-analysis 
methods  have  been  developed.  Some  versions  are  hybrids  belonging  to  more  than  one  of 
the  aforementioned  techniques. 

According  to  McPherson  (1976),  objective  analysis  methods  used  in  operational 
meteorology  can  be  divided  into  two  basic  categories.  The  first  represents  an  analyzed 
field  as  a  series  expansion  (spectral  analysis): 


Zg=aifi  +  a^2'*’ 


+  an/m 
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where  Zg  represents  the  departure  of  the  field  from  its  mean  value.  The  fj  represents  a  set 

of  orthogonal  functions—for  example,  a  cosine  series.  The  analysis  procedure  involves 
determining  the  time-dependent  coefficients  aj  which  make  the  series  expansion  best  fit 

the  observed  data  by,  for  example,  a  least-squares  technique.  Analysis  methods  of  this 
type  are  in  operational  use  at  the  National  Meteorological  Center  (NMC)  (Flattery  1971; 
Hayden  1976)  and  the  United  Kingdom  Meteorological  Office  (UKMO)  (Dixon  1976), 
The  second  basic  method,  called  the  "gridpoint"  method,  includes  statistical 
interpolation  and  successive  correction  techniques.  Here,  the  analyzed  value  Z|  at  a 
discrete  point  g  is  given  by  a  linear  combination  of  observations  that  are  nearby  in  time 
and  space: 


Zg=aiZ°+a^°+ . +  a^;  2.2 

where  Z®  represents  observed  values  at  the  several  stations  within  some  predetermined 
influence  radius  of  point  g,  and  the  coefficients  aj  determine  the  influence  of  each 

observation  on  the  analyzed  value.  In  this  case,  the  analysis  involves  determining  the 
coefficients  aj  in  the  linear  combination  for  each  point  of  the  analysis  grid. 

Two  forms  of  representation  are  associated  with  these  two  basic  methods  of 
analysis:  the  discrete  form  and  the  spectral  form.  In  the  former,  the  analysis  is  a  set  of 
values  at  discrete  points  in  space  and  time;  in  the  latter,  the  analysis  is  represented  by  a 
series  expansion  such  as  Eq.  2.1.  The  spectral  form  has  been  used  primarily  for  global 
and  hemispheric  applications.  The  discrete  form  has  been  extensively  used  for  limited- 
area  applications  in  mesoscale  meteorology. 


15 


2.4  Initialization  Techniques 

The  objective-analysis  methods  described  in  the  previous  section  generally  do  not 
provide  balanced  mass  and  wind  fields  to  initiate  a  forecast.  Uncompensated  errors  in 
wind  and  pressure-temperature  observations,  interpolation  of  observations  to  model 
grids,  and  the  numerical  model's  inability  to  exactly  describe  the  atmosphere  are  the 
primary  sources  of  this  dynamical  inconsistency.  The  dynamical  imbalances  in  the  initial 
data  lead  to  the  generation  of  spurious  inertia-gravity-wave  oscillations,  or 
"meteorological  noise."  Primitive  equation  models,  unlike  geostrophic  models,  admit 
these  higher-frequency  gravity-wave  solutions  that  can  have  amplitudes  much  greater 
than  their  counterpart  in  the  real  atmosphere.  These  gravity  wave  oscillations  can  obscure 
the  lower-frequency  Rossby-mode  component  of  the  model,  which  constitutes  the 
meteorological  signal.  Early  numerical  modelers  called  attention  to  the  need  to  eliminate 
the  spurious  high-frequency  oscillations,  which  can  compromise  the  forecast  procedure. 
First  of  all,  these  fast-moving  gravity  waves  require  short  computational  time  steps; 
second,  they  can  seriously  interfere  with  very  short-range  forecasts  (<12  h);  and  third, 
they  can  impair  vertical  motion  and,  hence,  precipitation  forecasts  (Daley  1981). 
Therefore,  a  long-standing  approach  has  been  to  eliminate  or  effectively  reduce  the 
amplitude  of  these  fast-moving  inertia-gravity  waves  at  imtial  forecast  time.  This  process 
is  known  as  model  initialization.  Chamey  (from  unpublished  letter  to  Phillip  Thompson, 
12  February  1947)  provided  insight  for  rectifying  this  initialization  problem  by 
suggesting  that  one  should  modify  the  initial  state  or  modify  the  governing  equations;  that 
is,  use  filtered  models. 

The  first  and  simplest  approach  was  to  exclude  any  possibility  of  the  high- 
frequency  oscillations  by  using  a  "filtered  system"  such  as  the  balance  equations  which 
simply  reduce  the  model  dynamics  to  the  quasi-geostrophic  response.  However,  this 


approach  severely  restricts  the  model  dynamics,  which  results  in  very  poor  forecasts 
beyond  24  hours.  The  so-called  primitive  equations  account  for  more  atmospheric 
dynamics  and  can  yield  much  better  forecasts,  so  they  are  generally  used.  However,  the 
primitive  equations  do  allow  the  amplification  of  fast-moving  gravity  waves  which 
requires  some  modification  of  the  initial  conditions  to  achieve  the  desired  dynamical 
balance. 

Over  the  years,  many  initialization  methods  have  been  developed.  A  summary  of 
the  more  widely  used  techniques  are  presented  here.  For  a  more  detailed  review,  the 
reader  is  referred  to  Daley  (1991).  In  static  initialization,  the  data  are  adjusted  at  a  single 
time  level  to  conform  to  some  dynamical  constraints.  That  is,  certain  time  derivatives  are 
identified  as  vanishing  in  order  to  eliminate  or  reduce  the  generation  of  inertia-gravity 
wave  noise.  In  the  conventional  static  initialization,  a  standard  practice  is  to  first  analyze 
the  geopotential  field  using  pressure-height  data  and  use  wind  observations  to  estimate 
the  gradient  of  the  geopotential  using  the  geostrophic  relation.  The  analyzed  geopotential 
fields  on  pressure  surfaces  are  then  used  in  the  mass  balance  equation  to  obtain  the 
streamfunction  of  the  nondivergent  wind,  from  which  the  rotational  wind  component  can 
be  computed.  A  major  limitation  in  using  the  balance  equation  to  determine  a  rotational 
wind  for  initializing  a  primitive  equation  model  is  that  the  lack  of  a  divergent  wind 
component  insures  the  presence  of  gravitational  modes  (Haltiner  and  Williams  1980) . 

The  most  common  approach  of  initialization  in  intermittent  FDDA  is  normal  mode 
initialization,  which  achieves  dynamical  balance  using  the  normal  modes  of  the  linearized 
dynamical  equations.  The  direct  use  of  normal  modes  was  introduced  by  Dickinson  and 
Williamson  (1972).  They  proposed  that  the  amplitudes  of  the  unwanted,  fast-moving 
modes  be  set  to  zero.  Their  method  was  effective  in  suppressing  the  spurious  noise  in 
linear  nKxlels,  but  failed  in  the  nonlinear  case.  A  nonlinear  normal  mode  scheme  was 
independently  developed  by  Machenhauer  (1977)  and  Baer  (1977).  In  nonlinear  normal 
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mode  initialization,  the  tendency  of  the  unwanted  modes  are  set  to  zero,  versus  setting  the 
amplitude  of  these  modes  to  zero.  The  solution  of  the  nonlinear  equation  requires  an 
iterative  process.  Unfortunately,  this  nonlinear  scheme,  without  the  inclusion  of  diabatic 
effects,  suppressed  the  meridional  circulation  in  the  tropics.  Puri  and  Bourke  (1982) 
used  the  idea  that  the  tropical  divergent  circulations  driven  by  convection  influence  mainly 
the  low  frequency  gravity  modes.  Therefore,  they  excluded  these  low  frequency  modes 
from  the  initialization  using  a  frequency  cutoff.  Wergen  (1982)  introduced  another 
method  in  which  average  diabatic  heating  is  obtained  by  integrating  the  model  for  a  few 
time  steps  prior  to  initialization.  This  model-produced  diabatic  heating  is  then  included  in 
the  nonlinear  forcing  in  the  iterative  process. 

In  global  models,  (Andersen  1977;  Daley  1979;  Temperton  and  Williamson  1981; 
for  example),  after  the  normal  modes  of  the  model  are  computed,  the  high-frequency 
inertia-gravity  waves  can  be  removed  by  projecting  the  inertial  wind  and  mass  fields  onto 
these  normal  modes.  However,  in  limited-area  models  it  is  not  possible  to  detine  the 
horizontal  structure  of  the  normal  modes. 

Bourke  and  McGregor  (1983)  introduced  a  method  of  initializing  a  limited-area 
model  without  explicitly  computing  horizontal  normal  modes.  In  this  technique,  termed 
vertical  mode  initialization,  the  free  modes  of  oscillation  of  the  prediction  model  are 
identified  by  linearizing  the  equations  about  a  basic  state  of  rest.  This  linearization 
permits  a  simple  decomposing  of  the  three-dimensional  eigenvalue  problem  into  a  series 
of  two  dimensional  problems.  The  vertical  decomposing  leads  to  a  number  of 
characteristic  vertical  modes,  one  corresponding  to  each  discrete  level  in  the  model. 
Balance  conditions  on  the  horizontal  structure  equations  are  then  derived  for  each  vertical 
mode.  Filtering  conditions,  in  which  the  tendencies  of  divergence  and  ageostrophic 
vorticity  are  set  to  zero,  are  applied  to  derive  linear  diagnostic  equations  for  the  mass  and 
divergence  fields.  In  Fig.  2.2,  from  Bourke  and  McGregor  (1983),  graphs  of  surface 


Figure  2.2.  Surface  pressures  (mb)  for  individual  grid  points  for  the  first  12  h  of 
prognosis  using  an  uninitialized  analysis  (solid  curves),  after  initializing  two  vertical 
inodes  (dot-dash),  and  after  initializing  two  vertical  modes  (dashes-using  Dirichlet 
boundary  conditions  and  dots-using  Neumann  boundary  conditions)  (Bourke  and 
McGregor  1983). 


pressure  before  and  after  initialization  using  the  Australian  regional  primitive  equation 
model  are  shown.  It  is  obvious  that  the  initialization  procedure  successfully  removed 
noise  from  the  integrations. 

Temperton  (1988)  devised  a  method  of  applying  Machenhauer's  criterion  without 
requiring  the  computation  of  the  coefficients  of  the  individual  modes.  Figure  2.3  (from 
Temperton  1988)  shows  graphs  of  500-mb  geopotential  before  and  after  the  application 
of  Temperton's  implicit  normal  mode  initialization  for  the  Canadian  finite-element  model. 
This  method  is  also  successful  in  removing  the  high-frequency  noise. 

Although  appropriate  for  intermittent  FDDA,  the  normal  mode  initialization  is  a 
distinct,  separate  step  from  the  objective  analysis  and  usually  leads  to  changes  in  model 
parameters.  As  a  result,  the  initialized  analysis  may  no  longer  fit  the  observations  as 
closely  as  desired.  An  alternative,  known  as  dynamic  initialization,  has  the  advantage  of 
avoiding  the  complications  of  computing  normal  modes.  In  this  method,  which  was 
introduced  by  Miyakoda  and  Moyer  (1968),  observations  are  inserted  (intermittently  or 
continuously)  over  a  period  of  time.  The  method  has  the  added  advantage  of  simplicity 
and  can  balance  physical  processes  as  well  as  the  mass  and  wind  fields.  In  dynamic 
initialization,  the  model  equations  were  integrated  forward  and  backward  through  time 
under  controls  that  encourage  time  derivatives  to  become  small  compared  to  spatial 
derivatives,  which  in  turn  selectively  dampened  the  high-frequency  components  of  the 
solutions.  Early  versions  of  this  technique  required  several  repetitions  of  the  integration 
cycle  to  successfully  reduce  the  gravity-wave  oscillations  and,  as  a  result,  these  schemes 
were  computationally  expensive.  Another  disadvantage  was  that  the  slow  modes  of  the 
model  were  dampened. 

More  recently,  researchers  including  Bratseth  (1982),  Sugi  (1986),  and  Satomura 
(1988)  have  developed  schemes  which  are  more  computationally  efficient  and  have  more 
selective  damping  properties.  Figure  2.4  (fiiom  Sugi  1986)  depicts  graphs  of  the  gravity- 
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Figure  2.3.  Time  trace  of  height  field.  Heavy  line;  no  initialization.  Light  line:  after 
two  iterations  of  the  implicit  nonlinear  normal  mode  initialization  scheme  (Temperton 
1988). 


VERTICAL  MODE  (j) 


Figi^  2.4.  Gravity-wave  activity  (before  and  after  mitialization)  for  five  vertical  modes 
of  a  baioclinic  model  (Sud  1986). 
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wave  activity  (before  and  after  initialization)  for  five  vertical  modes  of  a  baroclinic  model. 
The  gravity-wave  noise  is  dramatically  reduced,  particularly  for  the  modes  of  large 
equivalent  depths.  Dynamic  initialization,  as  well  as  the  Laplace  transform  and  bounded 
derivative  methods  (described  below),  is  well  suited  for  initializing  data  on  a  limited 
domain. 

Lynch  (1985a)  developed  an  effective  method  of  initialization  based  on  a  filtering 
scheme  which  uses  a  modified  inverse  Laplace  transform.  This  technique  is  equivalent  to 
the  nonlinear  normal  mode  initialization  method,  but  it  has  the  advantage  of  not  requiring 
a  transformation  of  the  model  equations  into  normal  mode  space.  Therefore,  the  Laplace 
transform  method  is  well  suited  for  initializing  limited-area  models  with  complex 
boundary  conditions.  Lynch  (1985b)  used  his  method  to  initialize  data  for  a  barotropic 
limited-area  model,  successfully  removing  high-frequency  gravity  wave  oscillations 
during  the  model  integration.  This  technique  has  also  been  applied  in  a  filtering 
integration  scheme  for  continuous  data  assimilation. 

Another  technique,  which  has  been  used  to  initialize  models  of  limited  domain,  is 
the  bounded  derivative  method.  Kreiss  (1979)  developed  the  methodology  of  controlling 
the  amplitudes  of  the  high-frequency  inertia-gravity  waves  by  requiring  the  derivatives  of 
the  model's  dependent  variables  with  respect  to  time  to  be  bounded,  i.e.,  of  order  unity, 
at  the  initial  time.  First,  the  equations  of  motion  are  nondimensionalized  so  that  certain 
terms  are  multiplied  by  a  small  parameter,  e.  Then,  if  the  first  derivatives  are  bounded, 
the  equations  can  only  be  satisfied  if  the  model  atmospheric  flow  is  geostrophic  and 
nondivergent  to  order  e.  If  the  second-order  time  derivatives  are  bounded,  the  resulting 
diagnostic  relationships  are  the  quasi-geostrophic  omega  equation  and  the  nonlinear 
balance  equation. 
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2.5  Various  Data-Assimilation  Methods 

An  assimilation  method  extensively  used  in  operational  meteorology  is  the 
analysis-forecast  cycle  technique,  commonly  referred  to  simply  as  intermittent  data 
assimilation  (see  Fig.  2.5).  This  process  consists  of  four  steps  which  are  repeated  at 
each  assimilation  cycle  (typically  every  3-12  h).  After  the  data  have  been  checked 
(quality  controlled),  a  static  three-dimensional  objective  analysis  (typically  successive 
correction  or  optimal  interpolation)  is  performed  using  observations  and  a  background 
field.  The  background  or  "first  guess"  is  usually  a  prior  model  forecast  valid  at  the 
analysis  time,  or  can  simply  be  climatology  or  a  combination  of  both.  Then,  the  analyzed 
fields  are  adjusted,  or  initialized,  to  conform  to  some  dynamical  constraint(s)  in  order  to 
reduce  or  eliminate  inertia-gravity-wave  noise.  The  final  step  consists  of  a  short-range 
numerical  forecast  to  obtain  first-guess  fields  for  the  analysis  at  the  next  assimilation 
cycle.  Thus,  the  new  estimates  (analyzed  fields)  are  clearly  based  on  the  past 
observations,  being  carried  forward  in  time  by  the  model  forecast,  and  on  the  current 
observations.  The  intermittent  updating  process  is  entirely  appropriate  as  long  as  most 
available  data  are  taken  at  the  same  time,  for  example,  at  synoptic  times.  This  technique 
is  currently  used  at  most  of  the  world's  major  operational  forecasting  centers,  including 
the  NMC  (DiMego  1988),  the  Norwegian  Meteorological  Institute  (Gronas  and  Midtbo 
1986),  and  the  European  Centre  for  Medium-Range  Weather  Forecasts  (ECMWF) 
(Hollingsworth  1986). 

The  extensive  operational  use  of  intomittent  data  assimilation  is  primarily  due  to  its 
computational  efficiency.  In  addition,  this  method  normally  includes  a  normal  mode 
initialization  scheme  that  produces  a  balanced  mass/wind  initial  state.  A  disadvantage  of 
this  method  is  that  it  is  not  totally  suited  for  asynoptic  data  types;  that  is,  it  can  not 
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assimilate  data  continuously.  However,  updates  on  the  order  of  every  2  to  3  h  can  be 
made,  allowing  some  asynoptic  data  into  the  assimilation. 

Two  other  mathematically  elegant  methods,  Kalman-Bucy  filtering  and  the  adjoint 
method  using  variational  techniques,  have  emerged  as  state-of-the-art  methodology  in 
FDDA.  The  latter  is  briefly  reviewed  first.  For  a  more  rigorous  and  complete  discussion 
of  the  general  theory  of  these  two  techniques,  the  reader  is  referred  to  Daley  (1991). 

Variational  assimilation,  based  on  the  calculus  of  variations,  involves  the 
incorporation  of  dynamical  constraints  in  a  variational  treatment  and  has  been  pursued  by 
Sasaki  (1969),  Stephens  (1970),  and  others  for  several  years.  In  variational  calculus, 
stationary  points  (extrema)  of  integral  expressions  known  as  functionals  are  determined. 
7  is  a  functional  of  the  function  q(t)  in  the  interval  (tm,  tn),  if  it  depends  on  all  the  values 
q(t)  for  tm  ^  t  ^  tn-  In  this  approach,  successive  analyses  are  mutually  adjusted  to 
effectively  increase  the  database  at  each  time  step  by  using  information  at  other  analysis 
times  through  the  forecast  equations.  The  objective  is  to  produce  initialized  values  of  q 
subject  to  certain  constraints  such  as  the  hydrostatic  relation,  the  continuity  equation,  the 
geostrophic  relation,  or  the  nonlinear  baltmce  equation.  The  approach  is  designed  to  keep 
the  initialized  fields  close  to  the  observations  while  satisfying  the  constraint  (Daley  1991). 

Lewis  (1972)  developed  a  variational  scheme  using  a  thermal  wind  relationship  and 
the  hydrostatic  equation  as  constraints.  A  more  recent  variational  approach,  known  as  the 
adjoint  method  (Lewis  and  Derber  1985),  uses  a  complete  dynanuc  model  as  a  strong 
constraint.  This  method  fits  a  model  to  observational  data  distributed  over  a  finite  period 
by  computing  derivatives  of  model  output.  An  iterative  method  minimizes  the  weighted 
squared  difference  between  the  original  analyses  at  several  times  and  the  coincident 
solutions  to  the  model  (constraint)  for  a  given  output  variable.  The  final  analyses  are 
constrained  to  satisfy  the  model  forecast  from  a  set  of  initial  conditions.  The  functional, 
7,  is  minimized  by  finding  the  gradient  of  7  with  respect  to  the  initial  conditions. 
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The  output  of  any  model  depends  upon  a  set  of  input  variables:  initial  conditions, 
boundary  conditions,  and  even  modeling  parameters  of  physical  processes.  Because  of 
this  dependence,  the  "adjoint"  can  be  used  to  determine  the  sensitivity  of  the  model  output 
to  any  input  parameter,  i.e.,  initial  conditions.  Neglecting  forecast  errors  over  the 
assimilation  period  (time  period  spanning  the  observations),  this  technique  can  produce 
the  optimal  initial  state  (end  of  the  assimilation  period)  such  that  the  model  best  fits  the 
observations  through  the  entire  period. 

The  alternate  "next-generation"  FDDA  method  is  called  Kalman-Bucy  (K-B) 
filtering.  The  K-B  method  can  be  thought  of  as  a  continuous  dynamic  FDDA  where  the 
weighting  factors  are  optimally  determined  by  explicitly  calculating  the  error  covariance 
of  the  analysis.  In  1960,  Kalman  developed  the  basic  theory  for  the  linear,  unbiased  data 
assimilation  scheme  known  as  the  Kalman  filter  (Kalman  1960). 

The  K-B  filter  (Kalman  and  Bucy  1961;  Ghil  et  al.  1979  and  1981),  which  is  the 
time-continuous  counterpart  to  Kalman's  original  scheme,  uses  the  forecast  model  itself 
to  predict  the  background  error  statistics.  Hie  goal  of  this  method  is  to  obtain  the  most 
accurate  analysis  value  for  all  time  during  the  assimilation  period  using  only  present  and 
past  observations.  The  weight  given  to  the  current  observations  is  inversely  proportional 
to  their  variance,  and  the  accuracy  of  the  analysis  is  the  sum  of  the  accuracies  of  the 
forecast,  based  on  the  past  observations,  and  of  the  current  observations.  The  K-B  filter 
minimizes  the  analysis  error  variance  not  only  at  every  time  step,  but  over  the  entire 
assimilation  interval  in  which  data  are  provided  (Ghil  1989).  Through  an  application  of 
Bayesian  ideas  in  a  dynamical  sense  (Kalman  1960;  Lorenc  1986),  the  filter  is  able  to 
extract  aU  useful  information  from  the  observational  increment/residual  at  each  time  step, 
thus,  allowing  observations  to  be  discarded  as  soon  as  they  are  assimilated.  As  a  result, 
the  method  is  sequential. 
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The  K-B  filter  and  the  adjoint  method  are  very  promising  assimilation  techniques. 
Both  methods  have  been  shown  to  produce  improved  assimilated  states  for  model 
forecast  integration  as  compared  to  those  from  the  operational  techniques.  However, 
these  "next-generation"  data  assimilation  methods  are  of  limited  practicality  at  the  present 
because  of  their  complexity  and  extensive  computational  requirements.  Therefore, 
operational  implementation  is  still  several  years  away  (Stauffer  and  Seaman  1990). 
However,  Lorenc  (1988)  has  pursued  approximations  to  the  adjoint  technique  in  a  quest 
to  make  it  operational  in  numerical  weather  prediction. 

Most  research  has  been  directed  toward  the  less  elegant  but  more  practical  method 
of  dynamic  assimilation,  in  which  the  numerical  prediction  model  serves  as  an  integrator 
of  observations  distributed  in  time  and  space.  In  this  approach,  which  has  become 
known  simply  as  "nudging"  or  Newtonian  relaxation,  the  model  integration  is  interrupted 
periodically  and  the  current  model  state  is  updated  with  the  new  observations.  During  the 
assimilation  cycle,  or  preforecast  integration  period,  the  model  variables  are  gradually 
driven,  or  nudged,  toward  the  observations  by  extra  forcing  terms  in  the  equations 
(Anthes  1974;  Kistler  1974;  Hoke  and  Anthes  1976;  Davies  and  Turner  1977).  As  a 
result,  the  model  fields  are  gradually  corrected  and  no  further  dynamic  balancing  through 
initialization  is  required.  The  general  form  of  the  predictive  equation  of  variable  S  is: 

^  =  F(S,x,t)  +  G  W  (So-S).  2.3 

All  of  the  model's  physical  forcing  terms  (Coriolis,  advection,  etc.)  are  represented  by  F, 
where  S  represents  a  model  dependent  variable,  Sq  represents  observations  of  S,  x  is  the 

independent  spatial  variable,  and  t  is  time.  The  second  term  on  the  right  is  the  nudging 
term,  where  G  is  the  nudging  constant  (generally  10*3  to  10*^)  and  W  represents  a  four- 
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dimensional  weighting  function.  The  data  to  be  nudged  can  be  either  derived  or 
measured,  analyzed  to  a  grid  for  assimilation  into  the  model,  or  inserted  as  individual 
observations. 

This  technique,  which  has  been  successful  in  bringing  the  data  and  the  model  in 
harmony  and  providing  a  relatively  noise-free  start  for  the  forecast,  has  been  widely  used 
on  the  global  scale  (Lyne  et  al.  1982;  Krishnamurti  et  al.  1988)  and  on  the  regional  scale 
in  limited-area  models  (Anthes  1974;  Hoke  and  Anthes  1977).  Of  late,  several 
researchers  have  developed  a  new  interest  in  this  technique  (Stauffer  et  al.  1985;  Bell 
1986;  Ramamurthy  and  Carr  1987;  Kao  and  Yamada  1988;  Wang  and  Warner  1988;  Kuo 
and  Guo  1989;  Stauffer  and  Seaman  1987, 1990).  It  is  currently  used  operationally  at 
the  UKMO  for  both  global  (Lyne  et  al.  1982)  and  regional  (Bell  1986)  data  assimilation. 

The  nudging  technique  mentioned  above  does  possess  a  few  desirable  attributes; 
these  were  summarized  by  N.  L.  Seaman  (from  his  lecture  notes  presented  at  the  1990 
Summer  Colloquium  on  Mesoscale  Data  Assimilation,  Boulder,  Colorado):  1)  The 
assimilating  model  is  complete,  so  irreversible  processes  are  included  without  difficulty. 
2)  Any  data  type  that  can  be  represented  as  a  tendency  of  a  prognostic  variable  can  be 
assimilated.  3)  Observation-nudging  can  easily  assimilate  asynoptic  and  single-level 
data.  4)  Analysis-nudging  requires  that  the  analyses  be  performed  only  once  prior  to 
model  integration,  i.e.,  it  is  economical.  5)  Nudging  does  not  require  a  separate 
balancing/initialization  step.  6)  Nudging  is  conceptually  and  computationally  simple.  On 
the  other  hand,  this  method  has  a  few  disadvantages:  1)  The  nudging  constant  is 
generally  assigned  in  an  application-dependent  semi-arbitrary  manner.  2)  Observation- 
nudging  is  based  on  "continuous  analysis"  during  model  integration  and  can  become 
computationally  expensive.  3)  Analysis  nudging  is  not  well-suited  for  asynoptic  data 
types.  4)  Use  of  accurate  data  in  observation-nudging  may  cause  assimilation  of  local  or 
unrepresentative  components,  e.g.,  microscale  observations  spread  over  a  large  area. 
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2.6  Comparisons  of  Operational  Assimilation  Techniques 

Two  papers  from  the  recent  literature  (Ramamurthy  and  Carr  1987;  Kuo  and  Guo 
1989)  compared  different  data-assimilation  methods.  Unlike  the  purpose  of  the  present 
discussion,  both  of  the  aforementioned  articles  advocated  Newtonian  nudging.  The  sole 
purpose  here  is  to  offer  a  short  review  of  the  available  comparative  studies. 

Ramamurthy  and  Carr  (1987)  studied  the  applicability  of  several  assimilation 
techniques  currently  being  employed  in  operational  models.  A  prime  scientific  objective 
was  to  determine  the  "best"  way  to  assimilate  asynoptic  observations  in  limited-area 
models.  A  sequence  of  ten  assimilation  experiments  was  conducted  using  different 
update  procedures.  In  each  experiment,  their  limited-area  model  was  initialized  with 
ECMWF  FGGE  level  ni-b  data,  and  then  12-h  assimilations  were  performed  using  level 
n-b  data  from  the  1979  Summer  Monsoon  Experiment  (SMONEX).  Forecasts  were  then 
made  from  these  assimilated  states.  The  first  experiment  served  as  the  control  run  since 
no  assimilation  was  performed.  Figure  2.6  is  a  schematic  of  the  overall  assimilation- 
forecast  strategy;  four  types  of  data  assimilation  were  compared:  1)  In  the  static 
assimilation,  the  model  was  updated  only  once  at  the  end  of  the  preforecast  (assimilation) 
period.  An  initialization  step  was  then  taken  to  suppress  the  noise  associated  with  the 
external  inertia-gravity  mode;  the  internal  modes  were  not  initialized.  2)  The  intermittent 
data-assimilation  experiment  was  similar  to  the  static  case,  except  the  model  was  updated 
twice  during  the  assimilation  period.  3)  Four  experiments  were  conducted  using 
continuous  indirect  data  assimilation  in  which  the  model  was  updated  whenever  new 
observations  became  available.  4)  Newtonian  relaxation  was  used  in  the  last  three 
experiments;  the  model  state  was  nudged  toward  analyses  produced  from  the 
observations.  The  reader  should  see  Ramamurthy  and  Carr  (1987)  for  a  detailed 
description  of  the  experiments. 
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Figure  2.6.  A  schematic  of  the  assimilation-forecast  cycle  for  various  experiments 
(RamamurAy  and  Carr  1987). 


Comparisons  among  these  differing  techniques  were  made  by  examining 
assimilated  states  (analyses  obtained  at  the  end  of  the  12-h  assimilation  period). 
Newtonian  nudging  produced  better  assimilated  states  than  did  the  continuous 
assimilations  via  indirect  insertion.  Also,  the  continuous  assimilation  experiments 
produced  noisy  assimilated  fields  due  to  insertion  shocks. 

From  each  of  the  assimilated  states,  24-h  forecasts  were  made  and  the  results  were 
compared  against  each  other  and  with  the  observations.  The  continuous  assimilation  (via 
indirect  insertion)  continued  to  suffer  from  the  ill  effects  of  the  insertion  shock. 
However,  the  forecasts  from  the  nudging  experiments  had  a  minimal  amount  of  noise. 

The  degree  of  spinup  was  examined  in  terms  of  the  development  of  precipitation. 
The  excessive  shocking  associated  with  continuous  insertion  was  detrimental  to  the 
spinup  process  and  consequently  to  the  rainfall  predictions.  The  rotational  nudging 
experiment  (only  the  rotational  component  of  the  wind  was  nudged)  produced  the  most 
accurate  rainfall  predictions.  The  lack  of  divergence-related  shock  in  this  approach  aided 
in  the  evolution  of  the  physical  processes. 

The  predicted  tracks  of  the  cyclone  in  the  SMONEX  case  study  were  compared. 
During  the  first  12  h  (assimilation  period),  the  cyclone's  movement  was  erratic  in  all  of 
the  experiments  due  to  the  rapid  mutual  adjustment  between  the  mass  and  wind  fields  that 
occurred  in  the  early  stages  of  integration.  In  the  static  and  intermittent  assimilation 
experiments  in  which  the  fields  were  replaced  completely  at  the  end  of  the  assimilation 
period,  this  problem  extended  into  the  forecast  stage  for  another  12  h.  Overall,  the 
nudging  experiments  predicted  more  accurate  tracks  than  the  other  experiments. 

In  another  comparative  study,  Kuo  and  Guo  (1989)  examined  different  data- 
assimilation  strategies.  They  conducted  a  series  of  observing-system  simulation 
experiments  (OSSEs)  to  test  a  Newtonian  nudging  technique  for  continuous  assimilation 
of  observations  fi'om  a  hypothetical  network  of  profilers.  Twenty  experiments  were 
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described  in  their  paper.  Similar  to  the  woric  of  Ramamurthy  and  Carr,  they  investigated 
static  initialization,  intermittent  assimilation,  and  Newtonian  nudging. 

Comparisons  of  the  three  assimilation  techniques  were  made  by  obtaining  time 
series  of  wind  and  temperature  errors  during  the  assimilation  period.  With  nudging,  the 
errors  decreased  gradually  during  the  period.  The  intermittent  assimilation  produced  a 
stepwise  decrease  in  the  errors  and  the  results  after  four  assimilation  cycles  were 
considerably  better  than  in  the  static  assimilation  case.  This  reveals  that  intermittent  data 
assimilation  was  also  effective  in  producing  an  improved  initial  state  for  the  model 
forecast. 

When  considering  model  noise  generated  during  the  preforecast  period,  a  strong 
difference  emerged  between  continuous  nudging  and  the  interrrtittent  assimilation.  In  the 
nudging  experiment,  the  model  noise  gradually  decreased  during  the  assimilation  period. 
In  contrast,  the  intermittent  data  assimilation  produced  model  noise  with  large  spikes 
immediately  after  each  reanalysis  of  the  wind  field.  However,  this  noise  could  have  been 
reduced  substantially  if  a  normal  mode  initialization  procedure  had  been  incorporated. 

The  two  papers  reviewed  above  suggest  that  Newtonian  nudging  is  the  best 
assimilation  method  among  those  compared.  However,  the  investigators’  findings  may 
be  scientifically  inconclusive  due  to  their  incomplete  tests  (notably,  their  failure  to  include 
an  effective  normal  mode  initialization  procedure  in  the  intermittent  data  assimilation 
experiments).  Furthermore,  their  comparisons  were  limited  to  operational  data 
assimilation  techniques.  The  next-generation  methods  will  produce  more  optimal 
assimilated  states  for  initiating  a  numerical  weather  prediction  model.  But  as  stated 
previously,  these  methods  are  still  in  the  developmental  research  stage  and 
computationally  efficient  versions  are  not  cuirently  available  for  operational  use. 
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2.7  FDDA  Research:  Implications  and  the  Future 

FDDA  has  been  developed  tremendously  during  the  past  15  years  and  is  now  an 
essential  component  of  numerical  analysis  and  prediction  systems  in  both  research  and 
operations.  The  improving  performance  of  medium-range  global  and  hemispheric 
prediction  and  the  rapid  development  of  limited-area  and  mesoscale  models  at  centers 
such  as  the  NMC,  ECMWF,  and  UKMO  clearly  illustrate  the  practical  benefits  of 
research  in  data  assimilation. 

The  systematic  development  of  data  assimilation  methods  has  made  possible  the  use 
of  unconventional  and  asynopdc  observations  from  sateUites,  aircraft,  drifting  buoys,  and 
soon  fipom  wind  profilers  and  doppler  radars.  As  a  result,  the  accuracy  of  short-range 
forecasts  has  signiHcantly  improved,  although  only  minor  clianges  have  occurred  in  the 
global  observing  system  since  the  Global  Weather  Experiment  in  1979  (Bengtsson  and 
Shukla  1988).  According  to  Lange  and  Hellsten  (1986),  the  3-day  rms  forecast  error  for 
the  Northern  Hemisphere  dropped  by  more  than  35%  between  1979  and  1986.  During 
the  same  period,  medium-range  weather  prediction  has  been  usefully  extended  in  the  time 
scale  from  3-4  days  to  about  7  days  in  the  Northern  Hemisphere  and  to  4  days  in  the 
Southern  Hemisphere  (Bengtsson  1985;  Bourke  et  al.  1985).  The  largest  improvement 
has  occurred  at  middle  and  high  latitudes  of  the  Northern  Hemisphere.  Weather 
prediction  in  the  tropics  has  not  improved  nearly  as  much,  due  to  insufficient 
observations  and  deficiencies  in  the  formulation  of  critical  physical  processes. 

Lorenz  (1982)  addressed  the  limit  of  medium-range  predictability  and  concluded 
that  it  is  possible  to  predict  instantaneous  weather  patterns  with  better  accuracy  than 
guesswork  nearly  two  weeks  in  advance.  Such  extensions  in  predictability  will  depend 
on  the  ability  to  improve  current  data-assimilation  systems  as  well  as  the  numerical 
prediction  models  themselves. 
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Research  is  ongoing  to  find  ways  to  improve  data  assimiL  don  methods.  A  number 
of  research  centers  are  investigating  higher-resolution  models,  better  numerical 
techniques,  improved  physical  parameterizations,  and  improved  models  of  auto-  and 
cross-correlation  functions  for  prediction  error.  Another  area  of  interest  is  in  specifying 
the  diabatic  heating  at  the  initial  forecast  time.  A  detailed  specification  is  necessary  in 
order  to  correctly  analyze  and  forecast  the  divergent  wind  field  over  the  tropics.  A  few 
universities  are  making  major  strides  in  developing  the  next-generation  FDDA  techniques 
and  enhancing  existing  methods.  For  example,  adjoint  methods  are  being  intensely 
investigated  at  the  University  of  Oklahoma;  K-B  filtering  at  McGill  University,  Montreal, 
Canada;  and  Newtonian  relaxation  at  the  Pennsylvania  State  University. 

Several  other  critical  research  topics  need  attention  during  the  1990s.  They  include: 
assimilation  of  moisture  processes  (rainfall  data  and  integrated  liquid  water),  lateral 
boundary  conditions  and  mesoscale  predictability  (predictive  skill  of  limited-area  models 
is  strongly  controlled  by  the  lateral  boundary  conditions),  assimilation  of  surface 
conditions/characteristics,  and  quality  control,  which  will  continue  to  be  a  priority  as  new 
data  types  are  applied  in  FDDA.  Data  assimilation  is  relatively  new  on  the  meso-alpha 
and  sub-alpha  scales  and  needs  considerable  research  with  improved  datasets.  Existing 
first-generation  assimilation  systems  must  be  improved  to  effectively  forecast  on  these 
small  scales. 

During  the  next  decade,  FDDA  will  play  a  role  in  solving  or  alleviating  important 
environmental  issues  such  as  air  pollution  and  climate  control  (increasing  CO2  and  the 

ozone  hole).  Recent  clean-air  legislation  requiring  the  use  of  low-sulfur  fuels  and 
extremely  expensive  equipment  will  make  accurate  numerical  models  of  atmospheric- 
chemistry  transport  and  removal  more  important  than  ever.  FDDA  can  be  used  to  obtain 
accurate  meteorological  fields  for  input  to  complex  air-chemistry  nnodels.  Bengtsson  and 
Shukla  (1988)  have  suggested  that  a  comprehensive  analysis  of  global  observations 
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based  on  an  FDD  A  system  with  a  realistic  physical  model  should  be  used  to  produce 
internally  consistent,  homogeneous  datasets  for  the  earth's  climate  system.  These  global 
observations  will  include  many  new  environmental  variables  such  as  external  forcing 
variables,  concentrations  of  radiatively  and  chemically  important  rare  species,  and  land 
surface  and  oceanic  variables.  New  remote-sensing  systems  such  as  EOS  (Earth 
Observing  System)  will  monitor  these  variables. 

With  these  environmental  issues  in  mind,  Daley  (1991)  elaborated  on  a  possible 

vision  of  the  future  of  data  assimilation: 

"Firstly,  data  assimilation  will  no  longer  be  entirely  or  even  primarily 
concerned  with  short  (or  medium)  range  weather  forecasting. 
Secondly,  the  data  base  will  become  incredibly  diverse  in  both 
variables  measured  and  type  of  observing  system.  Thirdly, 
assimilating  models  will  be  much  more  comprehensive,  involving 
ocean,  land  surface  and  stratospheric  components.  Finally,  there  will 
be  considerably  more  emphasis  on  the  long  term  stability  of 
algorithms  to  facilitate  climate  change  signal  detection." 
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3.  NRL/NCSU  MESOSCALE  MODEL 

The  forecast  model  used  throughout  this  data-assimilation  research  study  was 
developed  at  the  Naval  Research  Laboratory  and  North  Carolina  State  University  and  is 
described  in  detail  by  Madala  et  al.  (1987),  Holt  et  al.  (1990),  and  Sashegyi  and  Madala 
(1990, 1992b).  The  use  of  the  full  model  is  based  in  the  philosophy  that  no  compromise 
should  be  made  in  any  components  of  the  assimilation  system.  The  model  generates 
forecast  of  temperature,  wind,  surface  pressure,  humidity,  geopotential  height,  vorticity, 
vertical  velocity,  and  precipitation. 

3.1  Dynamics  of  the  Model 

The  prediction  model  is  based  on  the  primitive  equations.  The  horizontal 
momentum  equations,  thermodynamic  energy  equation,  moisture  conservation  equation, 
and  pressure  tendency  equation  make  up  the  set  of  prognostic  equations.  The  hydrostatic 
and  continuity  equations  are  used  as  diagnostic  equations  for  geopotential  height  and 
vertical  motion.  The  prognostic  equations  include  forcing  terms,  diabatic  heating  terms, 
and  sink/source  terms,  which  are  regarded  as  physical  processes  and  are  treated 
separately.  The  primitive  equations  can  be  written  in  flux  form  in  the  sigma  coordinates 
as 
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with  the  hydrostatic  and  continuity  equations  given  by 
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where  the  horizontal  components  of  the  wind  u  and  v,  the  temperature  T,  specific 
humidity  q,  and  the  surface  pressure  p^  are  the  model  prognostic  variables.  The 

geopotential,  potential  temperature,  and  vertical  motion  in  sigma  coordinates  are 
represented  by  <{>,  0,  and  a,  respectively.  The  vertical  sigma  coordinate  (a)  is  defined 

by  the  ratio  of  the  pressure  p  to  the  surface  pressure  ps  (Phillips  1957),  hx  and  hy  are  the 
map  factors  for  the  horizontal  x  and  y  coordinates,  and  Pq  is  a  standard  reference  level  of 

1000  mb.  K  is  given  by  the  ratio  of  R,  the  gas  constant  for  dry  air,  to  Cp,  the  specific  heat 
at  constant  pressure  for  dry  air.  The  mass  divergence  D  is  defined  by 


V.(p,v)  =  J 


h,hy 


^(p,uhy)  +  ^(p,vhj 
ox  dy 


3.8 


and  D  is  the  vertically  integrated  mass  divergence  given  by 


3.9 


The  terms  Fu,  Fy,  Fj,  and  Fq  on  the  right-hand  side  (RHS)  are  forcing  terms  and 
represent  physical  processes  such  as  friction,  diffusion,  heating  and  moistening  due  to 
cumulus  convection  and  large-scale  precipitation,  among  others.  The  boundary  condition 
at  the  model  top  (a  =  0)  and  bottom  (a  =  1)  is  a  =  0. 

The  model  dynamical  equations  in  flux  form  can  then  be  written  in  finite  difference 
and  vector  notation.  Details  of  the  vector  and  matrix  elements  can  be  found  in  the  report 
byMadalaetal.(1987). 
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3.2  Numerics  of  the  Model 

This  hydrostatic,  primitive  equation  model  uses  a  terrain-following  coordinate 
system  and  has  a  one-way  interacting  nested  grid  network.  The  domains  of  the  fine  and 
coarse  grids,  with  the  smoothed  model  topography  used  on  each,  are  shown  in  Fig.  3.1. 
The  coarse  grid  domain  includes  the  continental  United  States  and  extends  from  40*W  to 
140*W  and  10*N  to  70*N  with  a  horizontal  resolution  of  2‘  longitude  (170  km  at  40*N) 
by  1.5*  latitude  (166.5  km).  The  GALE  (inner  or  fine)  grid  covers  a  smaller  domain 
including  the  eastern  half  of  the  United  States  and  extending  out  over  the  Gulf  Stream 
from  58*W  to  102*W  and  23.5*N  to  56.5*N  with  a  finer  horizontal  resolution  of  2/3* 
longitude  (56.7  km  at  40*N)  by  0.5*  latitude  (55.5  km).  With  this  grid  resolution,  the 
integration  time  step  used  is  180  seconds  on  the  coarse  grid  and  60  seconds  on  the  fine 
grid.  The  continuous  governing  equations  are  written  in  flux  form  and  the  finite 
difference  scheme  is  a  second-order  quadratic  conserving  scheme.  The  Arakawa  C-grid 
(Arakawa  and  Lamb  1977)  in  spherical  coordinates  is  used  for  horizontal  differencing. 
Lateral  boundary  conditions  suggested  by  Davies  (1976,  1983)  are  emphasized  in  the 
present  version  of  the  model  at  each  grid  boundary. 

In  the  vertical,  both  grids  use  10  equally  spaced  sigma  levels.  The  vertical  layers 
and  the  horizontal  staggered  C-grid  are  illustrated  in  Fig.  3.2  for  the  simple  case  of  only 
two  model  levels.  In  the  vertical,  the  layers  are  of  varying  thickness  Aok  with  the 
prognostic  variables  u,  v,  T,  and  q,  and  the  geopotential  ^  defined  at  sigma  levels  Ok- 
The  vertical  motion  a  in  the  sigma  coordinate  is  staggered  in  the  vertical,  being  defined  at 
sigma  levels  Ck+Aak/2  at  the  boundaries  between  the  vertical  layers. 


Figure  3.2.  Horizontal  staggered  C-grid  and  vertical  grid  network  utilized  in  the  model 
for  the  simple  case  of  two  model  sigma  levels  (K  =  2). 


3.2.1  Split-Explicit  Time  Integration 


The  model's  time  integration  scheme  is  the  split-explicit  method  which  allows 
larger  time  steps  by  effectively  separating  various  terms  in  the  prognostic  equations  into 
parts  governing  the  slow-moving  Rossby  modes  and  fast-moving  gravity  modes.  In 
applying  the  split-explicit  scheme  of  Madala  (1981)  to  the  nested  model,  the  time  step  on 
each  grid  satisfies  the  non-linear  CTL  condition  with  centered  differencing  for  the  phase 
speed  of  the  fourth  internal  vertical  gravity  mode  of  the  model.  Initial  estimates  of  the 
tendencies  of  the  model  variables  for  all  the  terms  are  computed  with  a  large  time  step  At 
These  explicit  first  estimates  of  the  tendencies  of  the  mass  (surface  pressure  p^  and 

temperature  T)  and  wind  components  (u  and  v)  are  then  corrected  for  the  motion  of  the 
higher  frequency  gravity  modes  (external  and  first  two  internal  modes). 

These  corrections  to  the  mass  and  momentum  tendencies  are  computed  by 
integrating,  at  shorter  time  steps,  the  amplitudes  of  the  deviations  of  the  divergence  and 
generalized  geopotential  from  their  values  at  time  t,  for  each  of  the  first  three  vertical 
modes.  Time  steps  of  At/8,  At/4,  and  At/2  are  used  for  the  external  modes  and  first  two 
internal  vertical  modes,  respectively.  The  average  of  these  deviations,  over  the  interval  of 
twice  the  large  time  step,  is  then  used  to  correct  the  initial  estimate  of  the  variables,  with 
the  non-linear  forcing  held  constant.  Madala's  split-explicit  scheme  is  extremely  cost- 
effective,  requiring  only  20-25%  of  the  computer  time  expended  by  explicit  methods. 
Plus,  Sashegyi  and  Madala  (1990)  showed  that  this  technique  reduces  the  amplitude  of 
the  freely-propagating  higher-frequency  gravity  waves. 
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3.2.2  Lateral  Boundary  Conditions 

To  update  the  flow  at  the  horizontal  boundaries  of  each  grid  during  the  integration 
of  the  NRL/NCSU  mesoscale  model,  boundary  conditions  are  required  for  the  horizontal 
components  of  the  wind  field  and  for  the  mass  and  humidity  fields.  Surface  pressure  ps, 
temperature  T,  and  specific  humidity  q  are  defined  at  the  lateral  boundary.  The  wind 
components  are  staggered  in  the  C  grid,  which  allows  the  tangential  wind  to  be  defined  at 
the  boundary  and  the  normal  wind  to  be  staggered  half  a  grid  point  in  fix)m  the  boundary. 
Idealized  boundary  values  and  tendencies  for  the  coarse  grid  are  derived  from  12-hourly 
NMC  2.5*  hemispheric  analyses.  The  fine-grid  boundary  values  are  obtained  from  the 
coarse  grid  integration. 

Computed  model  variables  are  relaxed  to  the  boundary  values  in  a  boundary  zone 
of  six  points  using  the  Davies  scheme  (1976, 1983).  The  model  variables  are  adjusted  at 
each  time  step  according  to: 

a  =  (l-a)a„,+  aab  3.10 

where  'a'  is  the  independent  model  variable  psu,  psv,  ps,  PsT,  or  psq,  and  is  a  linear 
combination  of  the  model  computed  value  am  and  the  prescribed  boundary  value  a^. 
Following  Gronas  et  al.  (1987),  a  is  defined  as  a  quadratic  function  of  the  distance  (n) 
from  the  lateral  boundary  in  units  of  the  grid  spacing: 

1  for  n  =  0 

6.5 -nl  ,  , 

a  =  — g —  for  l^n^5  3.11 

0  for  n>6 

Coarse-grid  boundary  values  are  computed  at  each  time  step  by  a  linear 
interpolation  in  time  from  the  12-hourly  NMC  analyses.  At  each  time  step,  model 
variables  on  the  coarse  grid  are  linearly  imopolated  to  the  six  points  in  the  boundary  zone 
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of  the  fine  mesh.  Then  these  values  are  linearly  interpolated  in  time  to  the  intermediate 
time  steps  for  the  fine  grid. 

3.3  Physics  of  the  Model 

Physical  parameterizations  suitable  for  simulating  extratropical  cyclones  are  used. 
These  model  physics  are  described  below. 

3.3.1  Planetary  Boundary  Layer 

The  planetary  boundary  layer  is  parameterized  using  a  single  layer  with  the  fluxes 
of  heat,  momentum,  and  humidity  computed  using  a  generalized  similarity  theory,  in 
which  the  drag  coefficients  are  stability  dependent  (Chang  1981).  The  exchange  of 
momentum,  moisture,  and  sensible  heat  between  the  surface  and  the  lowest  model  layer 
is  proportional  to  the  differences  in  the  wind  velocity,  humidity,  and  temperature  between 
those  two  levels  and  is  given  by: 


ap,v  _ 

ax 
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at 
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3.13 


where  x  is  the  stress,  H  is  the  sensible  heat  flux,  and  Xj  and  Hj  are  the  respective  surface 
values.  A  similar  relationship  holds  between  the  specific  humidity  q  and  the  moisture 
flux  E.  The  magnitudes  of  the  surface  stress  and  surface  fluxes  of  sensible  heat  and 
humidity  are  related  to  the  values  of  the  wind,  potential  temperature,  and  specific 
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humidity  in  the  lowest  model  layer.  In  this  parameterization,  only  positive  heat  and 
moisture  fluxes  are  computed.  Heat  and  moisture  loss  to  the  ocean  surface  is  not 
considered. 

3.3.2  Horizontal  Diffusion 

The  model  has  horizontal  diffusion  to  control  noise  accumulation  during 
integration.  The  diffusion  is  of  second  order  and  is  applied  directly  on  the  horizontal 
sigma  surfaces  using  a  forward  time  differencing  scheme.  The  model  variables  for  the 
flux  form  of  the  equations  are  smoothed  by  application  of  the  two  dimensional  diffusion 
operator  V»KhV.  The  operator  is  applied  directly  to  the  wind  components  PsU  and  psV, 
while  for  the  temperature  psT,  specific  humidity  psQ,  and  surface  pressure  ps,  the 
operator  is  applied  on  the  deviations  from  values  computed  at  the  same  sigma  levels  for 
the  U.S.  standard  atmosphere.  The  diffusion  coefficient  varies  with  latitude  as  the  square 
of  the  longitudinal  map  factor  hx,  so  that  the  coefficient  of  diffusion  decreases  as  the 
longitudinal  grid  distance  decreases  with  increasing  latitude.  We  define  the  diffusion 
coefficient  Kh  as 


where  p.  is  a  fixed  non-dimensional  diffusion  coefficient,  Ax  and  At  are  the  latitudinal 
grid  spacing  and  time  step,  respectively.  A  value  of  p.  =  0.005  was  found  to  provide 
sufficient  smoothing,  without  excessively  damping  the  synoptic  scale  waves  on  the 
coarse  grid.  On  the  fine  grid,  a  value  of  ^  =  0.006  is  used  to  provide  a  little  more 
damping  of  the  small  scale  features. 
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3.3.3  Cumulus  Convection 


Convective  heating,  precipitation,  and  moistening  of  the  environment  are 
parameterized  using  the  method  developed  by  Kuo  (1974)  and  modified  by  Anthes 
(1977).  Convective  heating  is  released  in  an  atmospheric  column  when  there  is 
conditionally  unstable  stratification  and  moisture  convergence.  Kuo  assumed  that  deep 
convection  occurs  in  a  conditionally  unstable  environment  where  there  is  large-scale 
convergence  at  low  levels  in  the  troposphere.  It  is  further  assumed  that  the  deep  clouds 
dissipate  and  mix  horizontally  with  the  environmental  air,  releasing  heat  and  moisture. 
The  total  moisture  used  to  form  the  clouds  in  a  given  area  is  given  by  the  total  horizontal 
convergence  of  moisture  into  that  area.  The  converged  moisture  is  partitioned  into 
moistening  and  heating  parts  where  the  ratio  is  determined  by  the  vertically  averaged 
relative  humidity  of  the  column.  The  vertical  heating  profile  is  determined  by  the 
difference  between  the  predicted  temperature  sounding  and  the  moist  adiabat  constructed 
from  the  boundary  layer  temperature  and  humidity. 

First,  the  horizontal  moisture  convergence  M  is  computed  for  all  model  layers 


using 
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The  total  horizontal  moisture  convergence  Mt  is  computed  by  summing  over  all  the 
layers  except  the  top  layer,  where  Mt  is  given  by 


3.16 


The  moisture  convergence  in  the  top  layer  is  computed  separately,  since  it  is  assumed 
that  clouds  do  not  penetrate  into  the  top  layer.  The  air  column  is  tested  for  convectively 
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unstable  lapse  rates  if  the  total  moisture  convergence  Mt  exceeds  a  critical  value.  The 
total  moisture  used  in  cloud  formation  is  dependent  on  the  vertically  averaged  relative 
humidity  rm  in  the  environment.  This  moisture  can  be  partitioned  into  a  fraction 
(1-b)  Mt,  which  is  condensed  and  precipitated  out  in  the  cloud,  and  a  fraction  b  Mt 
which  moistens  the  cloud  (or  condenses  and  re-evaporates  at  the  same  level),  where  the 
fraction  is  defined  as  b  =  l-rm-  Assuming  the  cloud  dissipates  and  mixes  with  the 
environmental  air,  then  the  moistening  and  heating  rates  of  the  environment  can  be  easily 
computed.  The  vertically  averaged  heating  rate  Qc  and  moistening  rate  Dc  satisfy 

(l_b)M.  =  317 

gLc 

bM.  =  .5£  3.18 

g 

where  L  is  the  latent  heat  of  condensation.  The  precipitation  Pc  from  cumulus  clouds 

V 

follows  from 

Pc  =  (l-b)gMt  3.19 

3.3.4  Large-Scale  Precipitation 

The  parameterization  of  large-scale  or  non-convective  precipitation/heating  fohows 
Manabe  et  al.  (1965);  that  is,  the  excess  moisture  in  a  supersaturated  layer  ^  ithin  a  stably 
stratified  environment  is  condensed  out  isobarically,  leaving  the  layer  saturated.  The 
evaporation  of  precipitation  in  unsaturated  layers  below  the  cloud  base  is  also  accounted 
for. 
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Following  the  outline  in  Haltiner  and  Williams  (1980),  if  q  and  T  are  the  initial 

values  of  the  specific  humidity  and  temperature  of  the  supersaturated  layer,  then  the 
resulting  saturation  specific  humidity  is  given  by 


q  +  Aq  =  qs(T+AT)  s  q,('I)  + 


BT 


AT 


3.20 


Jp 


where  Aq  and  AT  are  the  changes  (assumed  small)  in  the  specific  humidity  and 
temperature  of  the  layer.  Using  the  Clausius-Clapeyron  equation  to  define  the  saturation 
vapor  pressure  as  a  function  of  the  temperature  and  rearranging,  the  change  in  the  specific 
humidity  is  given  by: 


Aq 


q-qs 

l+Liqs/CpRvT2 
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where  Ry  is  the  gas  constant  for  water  vapor.  Aq  is  reduced  at  each  model  level  by  a 
factor  to  account  for  that  fraction  which  re-evaporates  in  the  next  layer  below.  The 
corresponding  AT  at  each  level  is  given  by  CpAT  =  -LcAq,  and  the  total  large-scale 
non-convective  precipitation  Pl  is  given  by: 


Pl= 


3.22 


3.3.5  Dry  Convective  Adjustment 

Also  following  Manabe  et  al.  (1965),  a  dry  convective  adjustment  is  applied  to 
remove  any  absolutely  unstable  lapse  rates  in  the  model  atmosphere.  In  the  real 
atmosphere,  unstable  lapse  rates,  which  can  develop  due  to  large  surface  fluxes  of 
sensible  heat  for  example,  are  removed  through  dry  convection  (plumes,  rolls,  etc).  In 
the  model,  these  mechanisms  are  sub-grid  scale.  However,  such  unstable  lapse  rates  will 
produce  unrealistically  large  vertical  velocities  in  the  model.  To  resolve  this  problem,  the 


vertical  column  at  each  grid  point  in  the  model  domain  is  tested  for  any  absolutely 
unstable  lapse  rates.  Then,  for  each  sequence  of  unstable  layers  in  the  vertical,  the  sub¬ 
grid  scale  vertical  mixing  is  parameterized  in  the  model  following  Manabe  et  al.  (1965). 
It  is  assumed  that  parcels  in  unstable  layers  are  lifted  adiabatically  to  the  next  layer, 
exchanging  their  heat  at  constant  pressure,  while  conserving  the  total  enthalpy.  The 
exchange  of  heat  takes  place  until  the  temperatures  of  the  unstable  layer  and  the  first 
stable  layer  above  the  unstable  layers  are  adjusted  to  a  neutral  lapse  rate. 

3.3.6  Surface  Parameters 

At  the  earth's  surface,  weekly-averaged  sea  surface  temperatures  (SST)  of  0.5' 
resolution,  taken  from  the  GALE  datasets,  are  interpolated  to  the  model's  fine  grid. 
Climatological  SSTs  for  the  months  of  January  and  February  of  one  degree  resolution, 
taken  from  Reynolds  (1982),  are  interpolated  temporally  to  obtain  a  weekly  mean  and 
interpolated  spatially  to  the  coarse  grid  of  the  model.  A  sea/ice  boundary  is  derived  from 
the  U.S.  Navy's  climatological  sea/ice  boundaries  for  month  of  January.  Model 
topography  for  both  grid  domains  is  obtained  from  Navy  10-minute  global  topographical 
data. 

3.4  Model  Initialization 

The  vertical  mode  initialization  scheme  of  Bourke  and  McGregor  (1983),  as 
applied  by  Sashegyi  and  Madala  (1990,  1992a)  and  Harms  et  al.  (1992a),  is  used  to 
reduce  the  undesirable  gravity  wave  oscillations  in  numerical  integrations  of  the 
NRL/NCSU  mesoscale  model.  The  scheme  applies  a  correction  to  the  Rossby  modes  by 
averaging  a  separate  integration  for  those  gravity  modes  whose  phase  speeds  are  much 
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larger  than  meteorological  systems.  In  this  scheme,  filtering  conditions  for  the  inertia- 
gravity  waves  are  applied  to  the  model  dynamical  equations  to  derive  linear  diagnostic 
equations  for  the  mass  and  divergence  Helds  for  the  first  few  vertical  modes  of  the 
limited-area  model.  The  process  sets  time  derivatives  of  the  gravity-wave  oscillations 
resulting  from  imbalances  between  the  initial  mass  and  momentum  fields  to  zero  by 
decomposing  the  total  model  tendencies  into  normal  modes.  With  suitable  boundary 
conditions,  the  equations  are  then  solved  iteratively  for  the  first  three  vertical  modes  of 
the  numerical  model  to  remove  the  high  frequency  gravity  waves  from  the  initial 
conditions.  The  scheme  has  been  shown  to  be  an  application  of  the  normal  mode 
initialization  scheme,  developed  by  Machenhauer  (1977)  and  Baer  (1977)  and  used  in 
global  models  (for  example  Andersen  1977;  Daley  1979;  Williamson  and  Temperton 
1981),  without  the  horizontal  structure  of  the  normal  modes  having  to  be  computed 
(Juvanon  du  Vachat  1986;  Temperton  1988). 

As  in  Bourke  and  McGregor  (1983),  the  mass  field  remains  fixed  at  the  lateral 
boundaries  of  the  model.  Bourke  and  McGregor  assumed  that  the  divergence  does  not 
change  at  the  boundary.  To  provide  a  more  realistic  boundary  condition,  an  approximate 
mass  divergence  is  computed  inside  the  lateral  boundary  using  the  thermodynamic 
equation.  In  the  present  scheme,  the  resulting  changes  to  the  tangential  and  normal  wind 
are  computed  assuming  no  change  in  the  stream  function  and  velocity  potential  along  the 
lateral  boundaries  of  the  model  domain.  By  using  the  C  grid  and  defining  the  tangential 
wind  at  the  lateral  boundary  and  the  normal  wind  staggered  half  a  grid  point  inside  the 
lateral  boundary,  the  changes  in  the  wind  field  along  the  boundary  are  then  naturally 
consistent  with  the  changes  in  the  integrated  vorticity  and  mass  divergence  over  the  model 
domain. 

Following  Temperton  and  Williamson  (1981),  the  vertical  structure  of  the 
model's  normal  modes  is  obtained  by  separating  the  linear  and  non-linear  terms  of  the 


governing  equations  and  then  linearizing  about  a  basic  state  at  rest  with  a  mean 
temperature  T*(ct),  which  varies  only  in  the  vertical.  A  generalized  geopotential  is 
defined  whose  derivative  gives  the  horizontal  pressure  gradient  in  the  horizontal 
equations  of  motion  and  whose  time  tendency  is  related  to  the  divergence:  A  separation 
of  the  vertical  structure  of  the  model  variables  is  then  possible.  The  model  dynamical 
equations  are  expressed  as  time  tendencies  of  mass- weighted  vorticity  mass  divergence 
D,  and  generalized  geopotential  O.  The  equations  with  the  P  term  included  with  the  non¬ 
linear  terms  on  the  RHS  are  written  as: 


^  +  fD  =  A^ 
dt  ^ 

3.23 

3D  2 

—  +  V  <D  -  fC  =  An 

3.24 

3t 

do 

—  +  X^D  =  A^ 

3.25 

at 

The  freely-propagating  inertia-gravity  waves  in  the  model  are  solutions  to  these  linearized 
equations  with  the  forcing  terms  on  the  RHS  of  the  equations  equal  to  zero.  Here,  f  is 
the  Coriolis  parameter,  D  is  the  mass  divergence,  u  and  v  are  the  horizontal  wind 
components,  and  O  is  the  generalized  geopotential,  defined  as: 


<D  = 


(})-(I)+RT-<l> 


3.26 


where  ps  is  the  surface  pressure,  <t>  is  the  geopotential,  (|)s  is  the  surface  geopotential,  R  is 
the  gas  constant  for  dry  air,  and  <|)*  is  the  average  geopotential  on  the  sigma  surface  of  the 
model.  Equations  3.23  -  3.25  describe  the  shallow  water  equations,  where  the  phase 
speed  of  the  free  gravity  modes  is  given  by  where  Xk  (which  can  be  represented  as 
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an  equivalent  depth  multiplied  by  the  acceleration  due  to  gravity)  is  the  eigenvalue  for  the 
k**'  mode. 

The  filtering  conditions  of  Bourke  and  McGregor, 


dt 


dt 


fC- V  <D 


=  0 


3.27 


on  the  divergence  D  and  the  ageostrophic  "vorticity"  f  ^  need  to  be  applied  to 

Eqs.  3.23  -  3.25  for  only  the  first  three  vertical  modes,  whose  phase  speeds  of  the 
gravity  modes  are  much  larger  than  about  25  m  s’^.  Since  the  terms  on  the  RHS  of  Eqs. 
3.23  -  3.25  depend  on  the  vorticity,  divergence,  and  geopotential,  the  resulting  set  of 
equations  is  solved  iteratively  (Bourke  and  McGregor  1983). 

The  resulting  changes  to  the  horizontal  wind  field  are  computed  from  the  changes  in 
the  mass  divergence  and  mass  weighted  vorticity  by  solving  the  following  Poisson's 
equations: 


2 


V  Ax  =  AD 

3.28 

2 

V  AY  =  a; 

3.29 

for  the  changes  in  the  velocity  potential  Ax  and  the  stream  function  A'P  for  each  sigma 

level  of  the  model.  Then  changes  to  the  mass  weighted  wind  follow  from: 

A(^u)  =  5^Ax-5yAy  3  3Q 

A(prv)  =  SyAx-f-M'I'  3.31 

Finally,  following  Temperton  (1984),  the  changes  in  the  surface  pressure  ps  and 


the  temperature  T  can  be  derived  directly  from  changes  in  the  generalized  geopotential. 
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4.  APPLICATION  OF  THE  BRATSETH  SCHEME 
FOR  THE  ANALYSIS  OF  GALE  DATA 

The  mesoscale  weather  prediction  model  described  in  Chapter  3  provides  a  shon- 
range  forecast  which  serves  as  the  background  field,  or  first  guess,  for  our  numerical 
analysis.  Following  a  brief  introduction  of  various  analysis  techniques,  an  in-depth 
description  of  the  NRL/NCSU  objective  analysis  scheme,  which  was  developed  and 
tested  during  the  course  of  this  research,  is  given. 


4.1  Introduction 

During  the  1960s  and  70s,  the  most  widely  used  objective  analysis  technique  in 
operational  data  assimilation  and  forecast  systems  was  the  method  of  successive 
correction.  This  type  of  analysis  was  first  introduced  by  Bergthorsson  and  Diios  (1955), 
and  a  number  of  versions  have  been  developed  by  Cressman  (1959),  Barnes  (1973),  and 
others.  It  is  still  widely  used  for  the  analysis  of  mesoscale  systems  in  the  research 
community  (Koch  et  al.  1983;  Benjamin  and  Seaman  1985;  for  example).  However, 
these  successive  correction  schemes  possess  a  couple  of  undesirable  attributes.  First,  the 
analysis  converges  to  the  data,  and  if  a  first  guess  is  used,  then  too  much  attention  is  paid 
to  the  observations  relative  to  the  first  guess.  Second,  observations  in  regions  of  high 
data  density  are  given  excessive  weight  relative  to  observations  in  areas  of  low  data 
density. 

In  the  1980s,  emphasis  shifted  in  operational  community  toward  the  more 
mathematically  complex  and  computationally  expensive  method  of  statistical 
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interpolation,  or  better  known  as  optimal  interpolation.  With  optimal  interpolation, 
analysis  errors  are  minimized  with  respect  to  the  spatial  structure  of  observational  and 
forecast  errors.  A  concise  review  of  objective  analysis  can  be  found  in  the  paper  by 
Schlatter  (1988).  Recently  the  method  of  optimum  interpolation  has  been  applied  with 
success  to  operational  limited-area  modelling  systems  (DiMego  1988;  Benjamin  1989; 
Mills  and  Seaman  1990). 

A  paper  by  Bratseth  (1986)  has  caused  a  recent  resurgence  of  the  successive 
correction  method.  He  introduced  a  successive  correction  scheme,  in  which  the  solution 
converges  toward  the  solution  obtained  by  optimal  interpolation.  His  technique  alleviates 
the  aforementioned  shortcomings  of  previous  successive  correction  methods  by  utilizing 
the  forecast  error  covariance  to  define  weights  which  are  reduced  in  regions  of  higher 
data  density  by  use  of  a  normalizing  factor.  Additionally,  the  Bratseth  method  requires 
much  less  computational  expense  than  does  optimal  interpolation  methods.  Therefore,  it 
is  very  attractive  for  use  for  a  limited-area  forecast  system  when  only  limited  resources 
are  available.  This  approach  to  successive  correction  methodology  is  currently  being 
used  operationally  in  a  multivariate  analysis  for  the  limited-area  forecast  system  at  the 
Norwegian  Meteorological  Institute  (Gronas  and  Midtbo  1986).  An  algorithm  similar  to 
the  Bratseth  scheme  has  been  applied  at  the  U.K.  Meteorological  Office  to  directly  insert 
data  into  both  the  global  and  regional  models  (Lorenc  et  al.  1991).  With  each  iteration  of 
the  scheme,  the  mass  and  wind  variables  are  updated  in  a  sequential  fashion.  To  maintain 
balance  in  their  models,  changes  to  the  geostrophic  component  of  the  wind  are  derived 
from  any  mass  changes;  the  wind  changes  are  non-divergent  and  the  divergence  during 
integration  is  damped. 

A  simplified  multivariate,  successive  correction  objective  analysis  scheme  using 
the  Bratseth  method  has  been  developed  at  the  U.S.  Naval  Research  Laboratory  (NRL) 
and  North  Carolina  State  University  (NCSU)  to  analyze  data  collected  during  the  Genesis 
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of  Atlantic  Lows  ExperimenL  During  intensive  observing  periods,  more  frequent  upper 
air  soundings  were  taken  at  selected  National  Weather  Service  (NWS)  sites  over  the 
region  of  the  eastern  U.S.,  supplemented  by  additional  sites  along  the  coast,  dropsondes 
deployed  from  aircraft  offshore,  and  a  denser  set  of  surface  observations  (Dirks  et  al. 
1988;  Raman  and  Riordan  1988).  As  has  been  the  practice  with  limited-area  forecast 
systems  using  optimal  interpolation,  the  length  scales  of  the  correlation  functions  in  the 
Bratseth  scheme  are  reduced  from  those  used  with  global  models  to  provide  more  detail  at 
the  smaller  scales  (DiMego  1988).  After  a  number  of  iterations  of  the  scheme,  the  length 
scale  of  the  correlation  functions  is  further  reduced  for  subsequent  iterations  to  speed  the 
convergence  of  the  scheme  (Gronas  and  Midtbo  1986;  Seaman  1988).  A  prior  forecast 
generated  by  the  NRIVNCSU  mesoscale  model  is  used,  as  in  the  operational  limited-area 
models,  to  provide  a  first  guess  for  the  analysis  scheme.  With  the  coarse  resolution  of 
the  upper-air  network,  it  is  only  by  using  such  a  mesoscale  model  to  generate  a  first 
guess,  is  there  any  hope  that  any  realistic  mesoscale  features  generated  by  the  model  will 
be  preserved  or  even  enhanced  by  the  analysis.  The  analysis  of  the  deviations  of  the  data 
from  the  first  guess  covers  a  region  similar  to  the  model's  fine  grid  but  at  a  lower 
resolution  of  1.5*  in  latitude  and  longitude.  The  grid  spacing  used  is  consistent  with  the 
average  upper  air  station  spacing  of  about  350  km  over  the  eastern  U.S.  To  remove  bad 
data  or  data  inconsistent  with  the  scales  resolved  by  the  numerical  model,  an  efficient 
quality  control  procedure  was  implemented  in  a  manner  similar  to  DiMego  (1988). 

In  this  chapter,  a  general  description  of  our  objective  analysis  method  is  given  and 
a  qualitative  evaluation  of  the  scheme  is  provided.  Experiments  with  the  analysis  scheme 
are  described.  The  vertical  mode  initialization  scheme  of  Bourke  and  McGregor  (1983)  is 
used  to  diagnose  the  divergent  circulations  associated  with  the  analyzed  mesoscale 
features  and  upper-level  jet  streaks.  The  interaction  of  these  secondary  circulations  in  the 
jet  streaks  can  be  important  for  cyclone  development  (see  Uccellini  et  al.  1984;  Uccellini 
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and  Kocin  1987;  for  example).  The  impact  of  these  mesoscale  circulations  produced  by 
the  analysis  on  subsequent  forecasts  is  tested  and  precipitation  forecasts  are  compared  to 
observed  values. 

4.2  The  NRL/NCSU  Objective  Analysis  Scheme 

A  simplified  application  of  the  Bratseth  (1986)  scheme  has  been  developed,  as 
compared  to  the  full  multivariate  implementation  used  at  the  Norwegian  Meteorological 
Institute  by  Gronas  and  Midtbo  (1986).  In  our  scheme,  univariate  analyses  of  the  mass 
and  wind  fields  are  first  performed  using  the  successive  correction  approach  of  Bratseth. 
The  deviations  from  our  first-guess  field,  which  is  derived  from  a  prior  12-h  forecast 
valid  at  the  analysis  time,  are  analyzed.  To  provide  the  coupling  between  the  wind  and 
mass  fields,  the  analyzed  wind  is  used  as  an  initial  estimate  of  the  geostrophic  wind  in  a 
further  iteration  on  the  geopotential.  The  geostrophic  wind  is  used  to  extrapolate  the 
geopotential  to  the  grid  point  locations  in  a  fashion  similar  to  Cressman  (1959),  but 
consistent  with  the  Bratseth  scheme.  In  each  subsequent  iteration,  an  improved 
geostrophic  wind  estimate  is  then  defined  by  the  new  geopotential  gradient,  and  not  by 
the  original  wind  as  is  widely  used  with  the  Cressman  scheme.  The  wind  analysis  is  then 
corrected  for  the  new  geostrophic  wind  as  in  Kistler  and  McPherson  (1975).  A  further 
analysis  of  the  wind  is  performed  to  improve  the  larger  scales  of  the  ageostrophic  wind. 
The  temperature  is  then  corrected  with  the  new  geopotential  thicknesses.  The  resulting 
analysis  corrections  are  interpolated  to  the  model  sigma  coordinates  to  update  the  first- 
guess  forecast,  so  as  to  not  destroy  any  balance  or  small  mesoscale  features  generated  by 
the  model. 

A  set  of  upper-air  soundings,  which  were  interpolated  to  10-mb  levels  from  the 
surface  up  to  1(X)  mb,  and  surface  data  for  the  second  Intensive  Observing  Period  (lOP) 
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of  GALE  were  obtained  fix)m  the  GALE  Data  Center  (GDC)  at  Drexel  University  (Mercer 
1987).  The  data  domain  included  the  eastern  U.S.  and  western  North  Atlantic  from 
115*W  to  45*W  in  longitude  and  from  10*N  to  60*N  in  latitude.  A  higher  density  of 
upper-air  soundings  and  surface  data  is  found  in  the  inner  GALE  region  covering  the 
Carolinas  and  southern  Virginia  from  the  Gulf  Stream  to  the  Appalachians.  An  analysis 
grid  covered  the  data  domain  with  a  resolution  of  1.5*  in  latitude  and  longitude  extending 
from  116*W  to  44*W  in  longitude  and  from  9*N  to  61.5*N  in  latitude.  To  correct  the 
model  forecast  first  guess  in  areas  with  no  data  over  the  ocean  and  around  the  lateral 
boundaries  of  the  analysis  grid,  bogus  data  derived  from  the  NMC  analysis  are  used. 
After  the  analysis  is  complete,  the  resulting  corrections  are  merged  with  corrections  to  the 
first  guess  derived  from  the  NMC  analysis  in  a  zone  along  the  lateral  boundaries  of  the 
analysis  grid.  The  corrections  are  then  bicubically  interpolated  to  the  model  fine  grid  and 
added  to  the  original  first  guess  on  the  pressure  surfaces  to  provide  a  final  analysis.  The 
corrections  themselves  are  also  interpolated  in  the  vertical  to  the  sigma  surfaces  of  the 
model  to  correct  the  original  first-guess  model  forecast.  By  applying  the  vertical  mode 
initialization  scheme,  both  the  divergent  and  the  ageostrophic  wind  are  derived  from  the 
analysis. 

A  further  important  part  of  any  analysis  scheme  of  real  data  is  the  quality  control 
of  the  data  prior  to  the  analysis.  Data,  which  are  inconsistent  with  the  first  guess  or  with 
neighboring  observations  or  are  of  a  scale  too  small  to  be  resolved  by  the  analysis 
scheme,  must  be  removed  to  prevent  aliasing  into  the  larger  scales  in  the  analysis. 

4.2.1  Data  Preparation  and  Quality  Control 

An  efficient  quality  control  procedure  has  been  developed  for  use  with  the 
analysis  scheme,  similar  to  those  at  operational  weather  centers.  The  soundings  are  first 
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checked  for  hydrostatic  consistency  and  smoothed,  retaining  data  at  50-mb  levels,  to  be 
more  representative  of  our  10-layer  model  vertical  resolution.  Similar  procedures  to 
those  used  at  NMC  in  the  regional  analysis  (DiMego  1988)  are  followed  by  performing  a 
"gross"  check  in  which  observations  with  large  deviations  fit>m  the  first-guess  values  and 
any  missing  data  are  removed.  In  a  simplified  "buddy"  check,  observations  not  agreeing 
with  neighboring  observations  are  removed.  Close  observations  are  replaced  by  an 
average  "super"  observation  and  any  remaining  isolated  observations  are  eliminated.  To 
smoothly  merge  with  the  NMC  analysis  at  that  synoptic  hour,  bogus  data  derived  from 
the  NMC  analysis  at  that  time  are  added  in  data  void  regions  such  as  over  the  ocean  or  at 
the  lateral  boundary  of  the  analysis  domain.  Over  the  ocean,  the  NMC  analysis  acts  as  a 
test  of  our  scheme.  A  further  advantage  of  the  quality  control  is  that  by  removing  bad 
and  missing  data  prior  to  the  analysis  our  analysis  code  is  able  to  be  vectorized  to  run 
rapidly  on  a  supercomputer.  Using  these  procedures,  the  influence  of  any  data  errors  and 
biases  is  reduced.  The  details  of  each  of  the  components  of  the  data  preparation  and 
quality  control  follow, 

a)  Hydrostatic  check  and  smoothing 

After  being  checked  for  hydrostatic  consistency,  the  soundings  are  smoothed  in 
the  vertical.  The  temperatures  and  the  u  and  v  wind  components  were  averaged  over  20- 
mb  intervals  weighted  by  the  log  of  pressure  for  the  temperatures  and  by  the  mass  for  the 
u  and  V  wind  components.  The  humidities,  which  were  more  noisy  (Kreitzberg, 
personal  communication),  were  smoothed  by  averaging  over  50-mb  intervals  using  the 
mixing  ratios  weighted  by  the  mass.  Consistent  geopotentials  were  then  recomputed  and 
variables  retained  only  for  levels  every  50  mb.  For  each  sounding,  a  "daily  average" 
surface  temperature  is  computed  by  linear  extrapolation  from  levels  above  the  boundary 
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layer.  The  sea  level  pressure  is  computed  by  extrapolation  of  the  average  surface 
temperature  to  sea  level,  and  a  plateau  correction  is  added  (Saucier  1955).  To  extrapolate 
heights  to  pressure  levels  below  ground  level,  an  average  of  the  surface  temperature  and 
the  daily  average  is  used.  For  efficiency,  the  soundings  are  sorted  into  5* 
latitude/longitude  boxes  for  each  pressure  level,  at  50-mb  increments  from  sea  level  to 
100  mb. 

b)  Gross  errors 

At  each  pressure  level,  data  are  checked  for  gross  errors,  as  compared  to  the  first- 
guess  analysis,  which  has  been  interpolated  to  the  pressure  surfaces  from  the  earlier  12-h 
model  forecast  In  this  gross  check,  those  upper-air  observations  that  differ  in  magnitude 
by  more  than  four  standard  deviations  fixrm  the  first  guess  are  rejected.  Any  missing  data 
specified  in  the  sounding  are  also  removed.  For  the  sea-level  pressure,  those  values  that 
differ  in  magtutude  by  more  than  two  and  a  half  times  the  standard  deviation  are  rejected. 
Estimates  of  the  standard  deviation  of  the  observations  from  the  first-guess  forecast  are 
obtained  by  adding  a  forecast  error  growth  to  the  observational  error.  The  values  used 
are  slight  modifications  to  those  used  at  operational  centers  (DiMego  1988;  Shaw  et  al. 
1987).  For  sea-level  pressure,  an  observational  error  of  1.5  mb  is  used  with  forecast 
error  growth  of  1.0  mb  in  6  h  to  give  a  first-guess  error  of  3.5  mb  for  a  12-h  forecast 
using  our  limited-area  model.  For  the  upper-air  variables,  the  values  used  for  the 
observational  error  and  the  forecast  error  growth  are  listed  in  Table  4.1.  For  the 
remaining  quality-control  procedures  and  the  analysis,  only  the  deviations  of  the  data 
from  the  first  guess  are  retained. 
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TABLE  4. 1 .  Standard  deviation  of  the  observational  errors  and  forecast  error  ^owth  for 
a  six-hour  forecast  period  in  parentheses.  These  are  the  input  error  statistics  for  the 
objective  analysis. 


Pressure 

Height 

Temperature 

Rel  humidity 

u,v  Wind 
(ms'^) 

(mb) 

(m) 

(*C) 

(%) 

100 

20  (17) 

2.0  (0.7) 

3.5  (1.3) 

150 

18(18) 

1.8  (0.6) 

4.0  (1.4) 

200 

15  (18) 

1.8  (0.5) 

4.0  (1.5) 

250 

14  (18) 

1.5  (0.6) 

4.0  (1.6) 

300 

14(18) 

1.0  (0.6) 

17(4) 

4.0  (1.6) 

350 

14(17) 

1.0  (0.5) 

15(4) 

4.0  (1.7) 

400 

12  (16) 

1.0  (0.5) 

13(4) 

3.5  (1.7) 

450 

10  (14) 

1.0  (0.4) 

13(4) 

3.5  (1.6) 

500 

9(12) 

1.0  (0.3) 

13(4) 

3.0  (1.6) 

550 

8(10) 

1.0  (0.3) 

13(3) 

3.0  (1.5) 

600 

7  (9) 

1.0  (0.3) 

12(3) 

3.0  (1.4) 

650 

6  (8) 

1.2  (0.3) 

12(3) 

3.0  (1.3) 

700 

6  (8) 

1.3  (0.3) 

12(2) 

2.5  (1.2) 

750 

6  (8) 

1.4  (0.3) 

12(2) 

2.5  (1.2) 

800 

6  (8) 

1.5  (0.4) 

13(2) 

2.5  (1.1) 

850 

6  (8) 

1.5  (0.5) 

13(2) 

2.5  (1.1) 

900 

5  (8) 

1.8  (0.5) 

13(2) 

2.5  (1.1) 

950 

5  (8) 

1.8  (0.6) 

13(2) 

2.5  (1.1) 

1000 

5  (8) 

1.8  (0.6) 

13(2) 

2.5  (1.1) 
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c)  Buddy  check 

A  simplified  buddy  check  procedure  is  used,  in  which  deviations  of  the 
observations  from  the  first  guess  are  compared  with  the  deviations  from  all  neighboring 
observations  including  itself.  Each  observation  is  compared  with  itself  and  all  the  others, 
within  a  maximum  distance  of  833  km  (equal  to  7.5*  of  latitude).  The  maximum  allowed 
deviations  for  a  pair  of  observations  is  a  linear  function  of  the  cartesian  distance  between 
the  observations,  and  is  equal  to  one  standard  deviation  at  zero  distance  and  increases  to 
3.5  times  the  standard  deviation  at  the  maximum  distance. 

d)  Average  close  observations  and  eliminate  isolated  observations 

To  speed  the  convergence  of  the  iterative  analysis  scheme,  observations  in  close 
geographical  proximity  are  combined  into  a  single  "super"  observation.  Each  5*  box  is 
divided  into  smaller  boxes  of  1.25*  longitude  by  1*  latitude,  and  an  average  of  all  the  data 
within  each  smaller  box  is  computed,  assigning  the  average  latitude  and  longitude  for  the 
position.  To  funher  prevent  isolated  data  from  adversely  affecting  the  analysis 
(generating  circular  corrections),  any  remaining  isolated  data  are  removed.  For  each  data 
box  with  less  than  two  of  the  surrounding  eight  boxes  having  any  data,  the  data  in  that 
box  are  removed. 

e)  Merge  with  NMC  analysis 

In  regions  of  no  data  over  the  ocean,  additional  "bogus"  soundings,  derived  from 
the  difference  between  the  first-guess  forecast  and  the  NMC  2.5*  hemispheric  analysis, 
are  used  to  merge  our  analysis  smoothly  with  the  NMC  analysis  over  the  ocean  at 
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synoptic  hours.  At  off-synoptic  times,  no  bogus  soundings  are  added  and  the  analysis 
returns  to  the  background  in  regions  of  no  data.  To  match  the  resolution  of  the  NMC 
analysis,  two  bogus  soundings  spaced  5/3*  apart  in  latitude  and  longitude  are  used  in 
each  empty  interior  box  requiring  bogus  data,  and  one  sounding  is  used  at  the  center  of 
each  boundary  box  requiring  bogus  data.  The  requirement  for  adding  bogus  data  to  an 
empty  box  depends  on  the  availability  of  data  in  the  neighboring  boxes.  In  this  case, 
additional  bogus  soundings  have  been  added  over  the  Gulf  of  Mexico  and  off  the  east 
coast  of  the  U.S.  An  example  for  1200  UTC  January  25  in  Fig.  4.1a  shows  the 
locations  of  the  geopotential  height  data  at  500  mb  after  the  quality  control.  Very  few 
upper-air  real  observations  (shown  by  the  crosses)  are  found  over  the  oceans,  where 
NMC  "bogus"  data  (shown  by  stars)  have  been  added.  For  the  sea-level  pressure, 
shown  in  Fig.  4.1b,  many  more  observations  were  available  over  the  ocean,  many  of 
which  are  superobs. 

4.2.2  Analysis  Interpolation 

I  he  attractive  features  of  the  Bratseth  successive  correction  scheme  (1986)  for  the 
interpolation  of  the  observations  to  the  analysis  grid  are:  1)  the  weights  are  dependent  on 
the  covariance  between  observations  aiid,  therefore,  are  dependent  on  the  accuracy  of  the 
observations,  2)  the  weights  are  reduced  in  regions  of  higher-data  density,  and  3)  the 
scheme  converges  to  that  obtained  by  optimal  interpolation.  As  has  been  done  in  the 
operational  implementation  of  this  scheme  in  Norway  by  Gronas  and  Midtbo  (1986),  the 
width  of  the  modelled  covariance  function  is  reduced  after  three  or  four  iterations  to  speed 
the  convergence  of  the  scheme  at  the  smaller  scales.  Further  passes  of  the  analysis 
scheme  on  a  smaller  grid  can  be  easily  added  in  regions  of  higher  observational  density, 
such  as  in  regions  of  drifting  buoys  or  dropwindsondes.  This  was  applied  to  the  case  of 


Figure  4.1.  The  location  of  the  data  for  (a)  the  geopotendal  height  at  SOO  mb  and  (b)  the 
sea  level  pressure  after  quality  control.  The  observations  are  indicated  by  the while 
bogus  values  are  indicated  by  The  analysis  grid  is  indicated  by  the  tick  marks  along 
the  domain  boundary. 


a  hurricane  with  an  earlier  version  of  the  analysis  based  on  the  Barnes  scheme  by  Shi  et 
al.  (1991).  The  steps  are  further  described  below. 

a)  Univariate  analysis  of  the  mass  and  wind  field 


In  the  first  step  of  the  analysis  scheme,  univariate  analyses  of  sea-level  pressure, 
geopotential,  temperature,  and  humidity  corrections  to  the  first  guess  are  successively 
adjusted  at  each  iteration.  In  the  Bratseth  scheme,  estimates  of  both  the  interpolated 
correction  and  the  observation  correction  are  computed  using  the  iterative  formulas.  For 
the  geopotential  corrections  <j)  to  the  first-guess  values  of  geopotential  provided  by  an 
earlier  12-h  forecast,  the  interpolated  value  at  the  grid  point  "x"  is  given  by: 


(v+l)  (v)  n 

<t>x  =  <>x 

j=l 

and  the  estimate  of  the  value  of  observation  correction  is  given  by: 


w 


4.1 


(tm-i)  |v)  n 
<l>i  =  <l>i  +  X«ij 
j=l 


'  '  tv) 

where  <t>j  is  the  value  of  the  observation,  and  (j)^  and  (|)j  are  the  interpolated  values  at 

the  grid  point  and  the  estimate  of  the  observation  value,  respectively,  for  the  \)th  iteration. 

(i)  (i) 

The  starting  corrections  <|)jj  and  <|)j  are  zero.  The  weights,  which  are  dependent  on  the 


o  (v)] 
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covariance  of  the  corrections  for  the  first  guess,  are  in  each  case  given  by: 
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where  the  covariances  arc  defined  in  terms  of  a  correlation  function  p  for  the  true  values 


2 

(j)t ,  and  is  the  ratio  of  the  observation  error  variance  Og  to  the  first-guess  (forecast) 
error  variance  a^.  For  the  forecast  error  variance,  we  use  the  same  variance  used  in  the 
quality  control  for  the  variance  of  the  observations  from  the  first-guess  forecast  listed  in 
Table  4.1.  A  "local  data  density"  mj  is  defined  at  the  observation  location  by 
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The  correlation  function  p(r)  for  the  geopotential  corrections  (|)  is  modeled  as  a  Gaussian, 


p(r)  =  4.6 

which  is  a  function  of  the  distance  r  between  observations.  The  correlation  length  scale 
"d"  is  defined  by  500  km  for  the  first  3  or  4  iterations  and  by  275  km  for  one  additional 
iteration.  Since  our  mesoscale  model  can  resolve  smaller  scale  mesoscale  features  than 
can  operational  regional  models,  a  somewhat  smaUer  length  scale  is  used,  as  compared  to 
that  used  at  NMC  (577  km)  by  DiMego  (1988).  After  the  first  pass  of  our  analysis 
scheme,  the  data  are  again  checked  for  inconsistency  with  the  first-pass  analysis  itself. 
Any  inconsistent  data,  which  differ  more  than  2.5  times  (1.5  times  for  sea-level  pressure) 
the  standard  deviation  from  the  first  pass,  are  removed  and  the  first  pass  is  repeated.  The 
same  correlation  function  is  used  for  the  univariate  analysis  of  sea-level  pressure, 
temperature,  and  relative  humidity.  With  observation  errors  included,  the  estimate  of 
observation  is  not  equal  to  the  analyzed  value  interpolated  to  the  observation  location,  but 
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it  converges  rapidly  to  the  value  of  the  observation  itself  (Bratseth  1986).  For  a  fixed 
length  scale  of  500  km,  the  convergence  of  the  observation  estimate  and  the  convergence 
of  the  analyzed  value  at  the  observation  location  are  shown  in  Fig.  4.2.  The  convergence 
is  measured  with  respect  to  the  number  of  analysis  iterations.  When  the  length  scale  is 
reduced  after  just  3  or  4  iterations,  the  resulting  rms  differences  for  both  the  observation 
estimate  and  analyzed  value  are  substantially  reduced,  as  depicted  in  Fig.  4.2.  It  can  be 
noted  however,  that  if  there  was  no  observation  error,  they  would  be  the  same. 

For  the  univariate  analysis  of  wind  corrections  to  the  first-guess  wind  field, 
different  correlation  functions  are  used.  The  correlation  functions  for  the  wind 
components  can  be  derived  from  the  correlation  function  p(r)  for  the  geopotential 
(Gronas  and  Midtbo  1986).  The  functional  form  of  the  correlation  functions  pu  and  pv 
for  the  components  of  the  wind  field 

Pu  = 

Pv  = 

where  p(r)  is  the  correlation  function  for  the  geopotential,  (x,y)  and  (xj.yj)  are  the 
cartesian  positions  of  two  observations.  These  correlations  effectively  give  smaller 
weights  across  the  direction  of  the  wind  than  along  it  without  having  to  compute  wind 
direction  as  in  Benjamin  and  Seaman  (1985).  A  modified  length  scale  "dy"  is  used  which 
is  different  from  the  length  scale  "d"  used  for  p(r).  A  value  of  700  km  is  assigned  for 
dy,  whereas  the  length  scale  for  p(r)  is  500  km.  This  value  of  dy  fits  the  length  scale 
found  by  Hollingsworth  and  Lonnberg  (1986),  where  the  transverse  correlations  in  the 
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Figure  4.2.  The  root  mean  square  (rms)  of  the  difference  between  the  observation 
estimate  and  the  observation  value  for  the  geopotential  at  500  mb  (lower  cuivc)  as  a 
function  of  the  iteration  number.  The  upper  curve  shows  the  rms  of  the  difference 
between  the  analysis  value  interpolated  to  the  observation  point  and  the  observation 
value.  The  nns  values  obtained  tdier  four  iterations  followed  by  a  further  iteration  with  a 
reduced  length  scale  are  indicated  by  tiie  X's. 
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wind  reach  zero.  After  3  to  4  iterations,  the  length  scales,  du  and  d,  are  reduced  to  380 
km  and  275  km,  respectively,  for  one  more  iteration  to  again  speed  convergence. 

b)  Geostrophic  enhancement  of  the  geopotential 


A  further  four  iterations  are  used  for  the  geopotential  correction,  using  improving 
estimates  of  the  geostrophic  wind  at  each  iteration  to  extrapolate  the  geopotential 
correction  to  the  grid  points.  An  initial  estimate  of  geostrophic  wind  is  obtained  from  the 

previous  univariate  u,  v  wind  corrections  by  linear  interpolation.  The  interpolated  value 

(v+i)  . 

<t>  at  the  grid  point  is  now  given  after  the  t)*"  iteration  by 


and  the  observation  estimate  (t)j  is 


+  f  a  '•’j  ~  ♦I  ■  -  yj)  {"i”  - 


where  <t>j  is  the  value  of  the  observation  and  is  an  estimate  of  the  geostrophic  wind 

th  W  ju.i) 

correction  at  the  observation  point  for  the  \)  iteration.  <j)j  and  Ug  j  are  the  pnor 

estimates  of  the  value  of  the  observation  of  geopotential  and  the  geostrophic  wind, 
respectively,  at  the  observation  point  For  the  next  iteration,  an  improved  estimate  of  the 

geostrophic  wind  correction  ^  **  obtained  at  the  observation  point  by  linear 
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interpolation  of  grid  point  values  computed  from  the  interpolated  values  of  the 


geopotential  at  the  grid  points. 


u 


(v+l) 

g.x 


-y  „ 
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u^^x  is  the  value  of  the  geostrophic  wind  at  the  grid  point,  derived  from  the  geopotential 

W  (u) 

<t)^  ,  and  Ug  j  is  the  value  at  the  observation  point  obtained  by  bilinear  inteipolation  from 
u^^x,  for  the  iteration,  "f  is  the  Coriolis  parameter,  and  f*  =  max(f(15*N),lf1) 
sign(f).  "c"  is  a  geostrophy  factor  depending  on  the  latitude  0  in  degrees  as 


1  for  101  >30* 

c  =  0.5  [  1  +  cos(120)  ]  for  15*  <  101  <  30*  4.13 

0  for  101  <  15* 


The  four  iterations  use  a  single  correlation  scale  "d"  of  500  km.  Only  the  larger  scales  are 
corrected,  without  changing  the  smaller  scales  which  are  not  as  geostrophic.  In  a  9-point 
zone  along  the  boundary  of  our  analysis  grid,  the  analyzed  values  are  merged  with  values 
interpolated  from  the  NMC  analysis.  The  analyzed  values  in  the  first  six  points  closest  to 
the  boundary  are  replaced  by  the  NMC  values,  while  a  linear  weighting  of  (.25,.75), 
(.5, .5),  (.75, .25)  is  used  to  merge  the  analyzed  values  with  the  NMC  values. 
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c)  Enhancements  to  the  wind  field  and  temperature  gradients 

Following  Kisder  and  McPherson  (1975),  the  geostrophic  pan  of  the  corrections  to 
the  wind  field  is  adjusted  to  match  the  new  corrections  for  the  geopotential.  That  is,  the 
univariate  wind  corrections  are  modified  by  adding  the  difference  in  the  geostrophic  wind 
between  the  geostrophic  and  univariate  analysis  of  the  geopotential.  Four  additional 
passes  of  the  univariate  wind  analysis  are  then  performed  to  enhance  the  corrections  of 
the  ageostrophic  wind.  A  similar  merging  with  wind  derived  from  the  NMC  analysis  is 
performed  in  the  9-point  boundary  zone  as  was  done  with  the  geopotential.  After  the 
corrections  have  been  computed  at  all  analysis  pressure  levels  of  the  analysis,  the 
univariate  temperature  corrections  are  adjusted  to  match  the  computed  thickness 
conections.  The  final  corrections  are  then  bicubically  interpolated  to  the  model  horizontal 
grid  and  added  to  the  first  guess  on  the  pressure  surfaces. 

The  analysis  deviations  are  interpolated  in  the  vertical  to  the  model  sigma  levels  and 
the  model  background  forecast  is  updated.  Before  interpolating  in  the  vertical,  an  initial 
estimate  of  new  surface  pressure  is  computed  from  the  analyzed  full  fields  on  the  p 
levels.  The  surface  pressure  of  the  background  forecast  is  similarly  computed  from  the 
fields  of  temperature  and  geopotential  height  on  the  pressure  levels.  The  difference  in 
surface  pressure  between  the  two  is  then  added  to  the  original  surface  pressure  of  the 
background  forecast  to  obtain  the  corrected  surface  pressure.  The  analysis  deviations  are 
then  linearly  interpolated  in  the  vertical  to  the  model  sigma  levels  and  added  to  the  first- 
guess  model  forecast.  The  thermodynamic  variables  (temperature,  humidity)  are 
interpolated  to  the  model  sigma  levels,  assuming  they  are  linear  in  the  log  of  pressure, 
while  the  wind  components  are  interpolated  assuming  they  are  linear  in  the  pressure. 
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4.3  Testing  the  Analysis  Scheme 

The  utility  of  an  analysis  scheme  may  be  judged  in  several  ways,  including  its 
ability  to:  1)  define  accurately  the  amplitude  and  phase  of  meteorological  waves,  2)  fit 
data  within  the  accuracy  of  observations,  3)  distinguish  between  accurate  and  erroneous 
data,  and  4)  provide  initial  conditions  for  numerical  weather  prediction  models  such  that 
forecast  error  is  minimized.  In  section  4.2,  it  was  shown  that  our  analysis  scheme  fits 
the  data  to  a  desirable  limit.  Here,  we  will  concentrate  on  items  1)  and  4)  from  a 
qualitative  viewpoint  by  comparing  our  analyses  and  subsequent  model  forecasts  based 
on  these  fields  with  the  mesoscale  analyses  of  Bosart  (1988)  and  Doyle  and  Warner 
(1990)  and  the  NMC  2.5“-resolution  hemispheric  analyses,  and  forecasts  produced  from 
these  NMC  fields  by  our  limited-area  model.  Application  of  the  vertical  mode 
initialization  scheme  on  both  the  NMC  and  our  analyses  is  used  to  compare  the  mesoscale 
vertical  motion  and  ageostrophic  wind  derived  fix)m  the  NRL  and  NMC  analyses.  Model 
forecast  precipitation  is  compared  in  each  case  to  observed  values.  Analyses  of  GALE 
raingauge  data  from  GDC  (Mercer  1987)  on  a  polar  stereographic  grid  of  about  84  km 
resolution  were  interpolated  to  a  0.5’  latitude/longitude  grid.  To  provide  values  outside 
of  the  region  of  the  GDC  precipitation  analyses,  precipitation  estimates  on  a  I’-resolution 
grid  derived  from  satellite  by  Martin  et  al.  (1988)  at  the  Cooperative  Institute  for 
Meteorological  Satellite  Studies  at  the  University  of  Wisconsin  were  also  interpolated  to 
the  0.5*  grid.  These  rainfall  estimates  were  derived  from  digital  infrared  images  of  the 
GOES -6  satellite  using  the  Global  Precipitation  Index  technique  of  Arkin  (1983).  This 
technique  assumes  that  rainfall  (R  in  mm)  is  linearly  related  to  the  fractional  cloud  cover 

(f): 


R  =  kxf 

where  the  coefficient  k  has  a  value  of  72  mm/day. 


4.14 
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4.3.1  Experimental  Design 

A  number  of  cases  were  performed  to  test  the  effectiveness  of  the  NRL/NCSU 
objective  analysis  scheme;  these  are  oudined  in  Table  4.2.  The  analyses  were  made  using 
surface  and  atmospheric  sounding  data  collected  during  lOP  2  of  GALE.  Twelve-hourly 
NMC  2.5*  hemispheric  analyses  for  the  period  0(X)0  UTC  25  January  to  1200  UTC  26 
January  1986  (25AX)Z  to  26/12Z)  were  used  for  the  initial  conditions  to  generate  the  first- 
guess  12-h  forecasts  and  for  comparison  to  our  analyses.  Six  model  variables  were 
analyzed:  sea  level  pressure,  geopotential  height,  u-  and  v-  wind  components, 
temperature,  and  relative  humidity.  Analyses  were  performed  at  19  vertical  levels  from 
1000  mb  to  100  mb  at  50-mb  increments  on  a  horizontal  grid  with  a  resolution  of  1.5* 
latitude  by  1.5*  longitude.  These  analyses  produced  from  the  NRL/NCSU  objective 
analysis  scheme  are  refetred  to  as  simply  NRL  analyses. 

In  case  4.1,  the  model  is  initialized  from  the  ISIWiL  NMC  analysis.  The  12-h 
forecast  serves  as  the  background,  or  first  guess,  for  univariate  (case  4.2)  and 
multivariate  (case  4.3)  analyses  using  the  NRL/NCSU  scheme.  For  each  case,  24-h 
forecasts  are  then  produced.  The  25/1 2Z  NMC  analysis  and  a  24-h  subsequent  forecast 
(case  4.4)  serve  as  a  baseline  for  comparison.  Cases  4.5  and  4.6  are  identical  to  cases 
4.3  and  4.4,  except  the  analyses  are  valid  at  26/OOZ  and  the  12-h  forecast  from  case  4.4 
serves  as  the  first  guess  for  the  NRL  analysis  in  case  4.5.  The  analyses  in  cases  4.7  and 
4.8  are  valid  12  h  later  at  26/12Z.  These  8  cases  essentially  serve  as  3  separate 
experiments:  1)  25/1 2Z  analyses  and  subsequent  forecasts,  2)  26/OOZ  analyses  and 
forecasts,  and  3)  those  at  26/1 2Z. 


TABLE  4.2.  Analysis  and  Forecast  Test  Cases 


Case 

4. 1  25AX)Z  NMC  analysis  and  12-h  model  integration 

4.2  25/12Z  univariate  analysis  with  12-h  forecast  from  case  4. 1  as  first 
guess  and  24-h  integration 

4.3  25/12Z  multivariate  analysis  with  12-h  forecast  from  case  4. 1  as  first 
guess  and  24-h  integration 

4.4  25/1 2Z  NMC  analysis  and  24-h  integration 

4.5  26/OOZ  multivariate  analysis  with  12-h  forecast  from  case  4.4  as  first 
guess  and  24-h  integration 

4.6  26/OOZ  NMC  analysis  and  24-h  integration 

4.7  26/1 2Z  multivariate  analysis  with  12-h  forecast  from  case  4.6  as  first 
guess  and  24-h  integration 

4.8  26/1 2Z  NMC  analysis  and  24-h  integration 
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4.3.2  The  Synoptic  Situation 

During  the  period  25/OOZ  to  27/1 2Z,  coastal  frontogenesis  and  cyclogenesis  took 
place  along  the  eastern  seaboard.  Riordan  (1990)  and  Doyle  and  Warner  (1990)  describe 
the  initial  stages  of  the  coastal  frontogenesis  from  24  January.  At  25/OOZ,  coastal 
frontogenesis  had  begun  offshore  along  the  Carolina  coast  in  an  onshore  easterly  flow 
regime  and  by  25/12Z  a  1000-km  long  coastal  front  paralleled  the  coastline.  Cold  air 
damming  east  of  the  Appalachians  and  warming  of  the  easterly  air  over  the  Gulf  Stream 
were  the  major  factors  contributing  to  the  strong  thermal  gradient  along  the  coast.  The 
easterly  flow  regime  can  be  seen  in  the  coarse  NMC  1000-mb  analysis  in  Fig.  4.3a, 
which  is  however  unable  to  resolve  the  strong  thermal  gradient  along  the  coast.  In  the 
upper  levels  of  the  atmosphere,  a  broad,  long-wave  trough  was  present  over  the  central 
U.  S.  (Fig.  4.3b).  An  associated  250-mb  jet  with  maximum  winds  of  55-60  m  s  * 
extended  from  northeast  Texas  across  the  Ohio  River  Valley.  During  the  next  24  h,  the 
trough  moved  eastward  and  deepened  rapidly  due  to  a  second  short  wave  digging 
southward  in  the  large-scale  trough.  By  26/0(^,  the  coastal  front  extended  from  Georgia 
to  southern  New  England  and  a  weak  low  moved  up  the  coast  along  the  front,  reaching 
the  Chesapeake  Bay  by  26/12Z  (Fig.  4.3c).  In  the  next  12  h,  the  coastal  low  merged 
with  a  weak  trough  seen  in  Fig.  4.3c  moving  eastwards  from  the  Great  Lakes. 

The  disturbance  which  ultimately  became  the  first  of  the  major  cyclones  of  GALE 
was  present  as  a  minor  frontal  wave  in  extreme  southeast  Texas  at  25/1 2Z.  The  surface 
disturbance  remained  weak  as  it  advanced  along  the  Gulf  coast  and  on  26  January  the 
surface  wave  weakened  further  as  it  moved  across  Georgia.  In  response  to  the  strong 
upper-level  short  wave  reaching  the  east  coast,  a  new  center  formed  off  Cape  Hatteras 
and  deepened  rapidly  as  it  moved  north  and  combined  with  the  pre-existing  waves 
associated  with  the  northward-moving  coastal  low  and  eastward  moving  trough.  By 


Figure  4.3.  (a)  Sea-level  pressure  (mb)  and  tempo-atuie  (’C)  and  winds  at  1000  mb,  and 
(b)  geopotential  height  (m).  temperature  (*C)  and  winds  at  500  mb  for  the  NMC  analysis 
on  2S/12Z.  (c)  and  (d)  as  in  (a)  except  for  26/12Z  and  27/12Z,  respectively.  Contours 
of  sea-level  pressure  and  temperature  at  1000  mb  are  4  mb  and  5*C,  respectively. 
Maximiun  wind  vectw  indicated  by  labelled  arrow  in  m  s'^. 
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27/1 2Z,  the  well-developed  low  pressure  center  crossed  the  New  England  coast,  as 
shown  in  Fig.  4.3d. 

4.4  Results 

Comparisons  between  first-guess  forecast  fields,  univariate  analyses,  multivariate 
analyses,  NMC  analyses,  and  Bosart’s  detailed  hand  analyses  are  made  at  25/1 2Z  and 
26/(X)Z.  Also,  model  forecasts  resulting  from  these  analyses  are  compared  with  each 
other  and  with  the  Bosart  analyses  valid  at  specific  forecast  times. 

4.4. 1  Comparison  of  25/1 2Z  Analyses 

Sea-level  pressure,  surface  temperature,  and  1000-mb  winds  for  the  first  guess 
(12-h  forecast  from  case  4.1),  and  the  NRL  analysis  at  25/12Z  are  shown  in  Figs.  4.4a 
and  b.  The  first  guess,  shown  in  Fig.  4.4a,  had  the  cold  air  damming  with  strong 
ridging  east  of  the  Appalachians.  The  presence  of  a  strong  temperature  gradient  along  the 
Carolina  coast  signified  the  existence  of  a  coastal  front.  The  NRL  analysis  shown  in  Fig. 
4.4b  retained  these  features,  suggesting  that  they  were  present  in  the  data.  The  25/1 2Z 
analysis  of  Doyle  and  Warner  (1990)  verifies  that  the  coastal  front  and  the  cold  air 
damming  were  present  (see  Fig  2b.  of  Doyle  and  Warner).  Note  that  the  coarse  NMC 
analysis  shown  in  Fig.  4.3a  is  unable  to  resolve  the  tight  thermal  gradient  along  the  coast. 
The  NRL  analyses  also  retained  the  strong  onshore  flow  from  the  first  guess  and  which 
is  readily  apparent  in  Doyle  and  Warner’s  analysis  along  the  Carolina  coast.  The  overall 
pressure  pattern  from  the  NRL  analysis  matched  Doyle  and  Warner's  analysis  much 
better  than  the  coarse  NMC  analysis  was  able  to  do.  Note  the  position  and  magnitude  of 
the  "low"  of  1000  mb  over  Lake  Superior  and  the  "high"  of  1046  mb  over  the  Gulf  of 


Figure  4.4.  Sea-level  pressure  (mb)  and  1000-mb  temperature  (*Q  and  winds  at  25/12Z 
for  (a)  the  first  guess,  a  12-h  forecast  from  case  4.1  and  (b)  NRL  analysis,  (c)  and  (d)  as 
in  (a)  and  (b)  but  for  the  geopotendal  height  (m),  temperature  (*Q,  and  winds  at  5(X)  mb. 
Contours  of  geopotendal  and  temperature  at  500  mb  are  60  m  and  2.5*C,  respecdvely. 
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St.  Lawrence.  The  first-guess  "low"  was  in  error  by  5  mb  and  was  corrected  totally  by 
the  NRL  analysis. 

The  analysis  of  500-mb  geopotential  height,  temperature  and  winds  are  compared 
to  the  first  guess  in  Figs.  4.4c  and  d.  The  NMC  500-mb  analysis  shown  in  Fig.  4.3b 
shows  a  cutoff  low  of  5627  m  over  the  western  Atlantic.  The  first-guess  field  provided 
by  the  model  forecast  failed  to  deepen  the  low  significantly,  and  only  a  weak  trough  is 
present.  The  univariate  NRL  analysis  (not  shown)  does  capture  this  feature,  but  the  low 
is  not  quite  deep  enough  (5634  m).  The  multivariate  analysis  (Fig.  4.4d)  is  able  to 
deepen  the  low  even  more  to  5625  m  by  extrapolating  geopotential  heights  using  wind 
information,  and  the  associated  winds  are  somewhat  stronger.  It  is  worth  noting  that  the 
data  used  in  the  NRL  analysis  in  the  vicinity  of  the  cutoff  low  are  in  fact  bogus  values 
derived  from  the  NMC  analysis  (see  Fig.  4.1a).  The  NRL  multivariate  analysis,  as  noted 
above,  reproduced  the  NMC  analysis  for  the  cutoff  low  with  remarkable  accuracy, 
considering  the  coarse  resolution  of  the  bogus  data  supplied  to  the  analysis.  The  first- 
guess  500-mb  temperature  field  is  out  of  phase  in  the  Great  Lakes  region.  The  -30'C 
isotherm  in  the  NMC  analysis  is  over  Lake  Michigan,  whereas  the  first-guess  forecast 
has  it  lagging  too  far  to  the  west  over  Minnesota.  Some  improvement  is  gained  with  the 
NRL  univariate  analysis,  but  the  phase  error  is  totally  corrected  by  the  multivariate 
analysis. 

The  NRL  analysis  is  successful  in  pushing  the  first-guess  500-mb  relative  humidity 

field  toward  the  NMC  analysis  (not  shown).  In  particular,  a  high  moisture  tongue  off  the 

Carolina  coast  was  absent  in  the  first-guess  field.  The  NRL  analysis  was  successful  in 

moistening  the  model  environment  to  the  appropriate  level.  The  250-mb  jet  associated 

with  the  upper-level  long  wave  trough  is  depicted  in  Fig.  4.5.  The  jet  from  the  first 

-1 

guess  is  too  weak  and  broad  with  a  jet  maximum  of  about  50  m  s  ,  compared  to  the 

-1 

NMC  analysis  which  has  a  jet  maximum  of  55  m  s  .  Both  the  univariate  and 


Figure  4.5.  Geopotential  height  (m),  temperature  (*C),  and  winds  at  250  mb  for  (a)  the 
first  guess,  a  12-h  forecast  from  case  4.1,  (b)  the  NRL  analysis,  and  (c)  the 
analysis  on  25/12Z  .  Contours  of  geopotential  are  every  120  m  and  isotachs  above  30 
m  S'*  are  shown  every  10  m  s"*. 
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multivariate  analyses  corrected  the  wind  speeds  and  sharpened  the  jet,  with  the 
multivariate  analysis  having  a  slightly  smaller  maximum  than  the  univariate  analysis. 
Acmally  the  analyzed  jet  maximum  in  the  multivariate  case  is  still  2  m  s  ^  stronger  than 
that  in  the  NMC  analysis.  However,  this  is  not  alarming  since  our  analysis  is  on  a  finer 
scale  than  the  NMC  analysis. 

The  first-guess  second-order  fields  of  absolute  vorticity  and  the  vertical  motion  are 
significantly  improved  by  the  NRL  objective  analysis.  The  500-mb  absolute  vonicity  and 
vertical  motion  fields  for  the  first  guess  and  NRL  analysis  are  shown  in  Fig.  4.6.  For 
comparison,  the  vorticity  and  vertical  motion  derived  from  the  NMC  analysis  are  shown 
in  Fig.  4.7.  The  vertical  motion  from  the  analyses  is  derived  by  application  of  the  vertical 
mode  initialization  procedure.  The  NRL  multivariate  analysis  greatly  improved  the 
overall  pattern,  including  the  gradients  of  vorticity,  when  compared  to  the  first-guess 
field.  The  vorticity  maxima  for  the  multivariate  analysis  are  close  in  value  to  those  in  the 
NMC  analysis.  In  the  first-guess  vertical  motion  field,  the  subsidence  west  of  the  cutoff 
low  is  strengthened  by  the  NRL  analysis  and  is  in  close  agreement  with  that  derived  for 
the  NMC  analysis.  However,  the  rising  motion  associated  with  the  developing  coastal 
front  in  the  first-guess  field  is  retained  in  the  NRL  analysis.  The  rising  motion, 
associated  with  the  surface  cyclone  located  on  the  southeast  Texas  coast  and  a  secondary 
circulation  of  the  upper-level  jet's  entrance  region,  is  evident  over  eastern  Texas  and 
northern  Mississippi.  These  two  centers  of  ascending  motion  are  better  defined  in  the 
NRL  analysis  than  either  the  first  guess  or  the  NMC  analysis. 

Cross  sections  were  produced  to  further  investigate  the  secondary  ageostrophic 
circulations,  associated  with  the  upper-level  jet  and  the  coastal  front.  The  locations  of  the 
cross  sections  are  displayed  in  Fig.  4.5c.  Cross  section  CSl  (94*W,  46*N  to  87.3*W, 
26*N),  which  extends  from  Minnesota  to  the  central  Gulf  of  Mexico,  cuts  across  the 
entrance  region  of  the  upper-level  polar  jet.  The  second  cross  section  (CS2  from  90*W, 
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Figure  4.6.  Absolute  voiticity  (in  units  of  s*0  with  geopotential  height  (m)  at  500 
mb  on  25/12Z  fcnr  (a)  the  first  guess  12-h  forecast  from  case  4.1,  and  (b)  the  NRL 
analysis.  The  vertical  velocity  motion  in  mb  h'^  at  S(X)  mb  for  (c)  the  first  guess  and  (d) 
the  NRL  analysis.  Contours  of  vorticity  ate  shown  for  6  units  and  above  every  2  units. 


Figure  4.7,  (a)  Absolute  vorticity  (in  units  of  lO'S  s'*)  with  geopotential  height  (m)  and 
(b)  the  vertical  velocity  in  mb  h'l  at  5(X)  mb  on  2S/12Z  for  the  NMC  analysis.  Vertical 
motion  of  magnitudes  -20  mb  h*l  or  greater  are  contoured  every  2,5  mb  h'^ 
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43.5*N  to  68.7*W,  26’N)  is  perpendicular  to  the  Carolina  coast.  In  all  figures  displaying 

cross  sections,  the  positions  of  the  grid  points  every  0.5*  in  latitude  are  denoted  by  the 

inner  tick  marks  for  facilitation  in  locating  meteorological  features  of  interest.  The  wind 

speed  normal  to  each  cross  section  is  shown  in  Fig.  4.8.  In  the  cross  section  across  the 

Gulf  of  Mexico,  the  NMC  analysis  at  25/1 2Z  (Fig.  4.9b)  shows  a  broad  region  of 

ascending  motion  in  the  right  rear  flank  of  the  polar  jet.  The  NRL  analysis  (Fig.  4.9a) 

has  enhanced  the  secondary  circulation  in  this  region.  The  maximum  vertical  velocity  is 

-1  -1 

40%  greater  than  in  the  NMC  analysis  (-14  mb  h  vs.  -10  mb  h  ).  Also,  more  evident 

in  the  NRL  analysis  is  the  ascent  associated  with  the  subtropical  jet  stream  to  the  south. 

Figure  4.10  depicts  the  cross  section  perpendicular  to  the  Carolina  coast  for  both  the 

NRL  and  NMC  analyses.  As  shown  in  Fig.  4.10a,  the  NRL  analysis  is  successful  in 

maintaining  the  shallow  ascent  associated  with  the  coastal  fiont,  which  was  generated  in 

-1 

the  first-guess  forecast.  Note  the  maximum  ascent  of  -1 8.3  mb  h  about  15(X)0  km  from 
the  northernmost  end  (left  side  of  cross  section)  in  Fig.  4,10a.  This  circulation  is  much 
stronger  and  more  well-defined  than  the  corresponding  circulation  firom  the  NMC 
analysis  (Fig.  4.10b).  The  subsidence  west  of  the  cutoff  low  in  the  subtropical  jet  stream 
and  ascent  east  of  the  approaching  upper-level  trough  (see  Fig.  4.6)  are  also  seen  s  both 
cases.  Increased  humidities  at  low  levels  offshore  and  over  the  coastal  front,  with  dryer 
values  behind  the  firont,  are  seen  in  Fig.  4.10c  for  the  NRL  analysis.  A  narrower  tongue 
of  high  humidity  is  also  found  at  mid-levels. 

4.4.2  Comparison  of  Model  Forecasts  (25/12Z-26/12Z) 

Forecasts  using  the  different  analyses  to  derive  the  initial  conditions  are  compared 
in  cases  4.3  and  4,4.  After  6  h  of  model  integration  (at  25/1 8Z),  the  forecast  based  on 
the  NMC  analysis  (case  4.4)  still  contained  a  relatively  weak  thermal  gradient  in  the 
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Figure  4.8.  Magnitude  of  the  wind  normal  to  the  plane  of  cross  sections  (a)  CS 1  and  (b) 
CS2  at  25/1 2Z  for  the  NRL  analysis  (see  Fig.  4.5c).  Contours  every  5  m  s*^  Potential 
temperature  (*K)  is  also  shown  by  the  dashed  contours  every  5*K  in  the  troposphere.  The 
position  of  the  grid  points  every  0.5*  in  latitude  are  shown  by  the  inner  tick  marks. 
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Figure  4,9.  The  vertical  motion  and  ageostix^hic  wind  in  the  plane  of  cross  ^tion  CS 1 
at  25/1 2Z  for  (a)  the  NRL  analysis  and  (b)  for  the  NMC  analysis.  The  magnitude  of  the 
vertical  motion  is  represented  by  the  solid  contours  at  -40,  -25,  -20,  -15,  -12.5,  and 
every  2  mb  h*l  for  values  -10  and  above. 


Figure  4.10.  The  vertical  motion  and  ageostrophic  wind  in  the  plane  of  cross  section 
CS2  at  25/12Z  for  (a)  the  NRL  analysis  and  (b)  for  the  NMC  analysis,  (c)  and  (d)  as  in 
(a)  and  (b)  but  for  the  relative  humidity,  contoured  every  10%  for  30%  and  above. 
Contours  of  vertical  motion  as  in  Fig.  4.9.  See  Fig.  4.5c  for  location  of  CS2. 
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Figure  4.10.  Continued. 
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vicinity  of  the  coastal  front  (shown  in  Fig.  4.11a),  since  the  model  was  still  in  a  spin-up 
mode.  An  inverted  surface  trough  is  present  near  Jacksonville,  Florida  in  Fig.  4. 1 1  b  in 
the  case  using  the  NRL  analysis  (case  4.3).  This  feature  is  weaker  and  further  east  in  the 
NMC  case.  Bosart's  analysis  (Fig.  4.12)  verifies  the  presence  of  the  trough  with  the 
genesis  of  a  cyclonic  circulation.  A  weak  surface  low  is  also  indicated  just  north  of  New 
Orleans  in  better  agreement  with  the  Bosart  analysis,  than  the  NMC  case.  The  weak 
coastal  circulation,  which  forms  a  low  pressure  center  that  moves  northeast  along  the 
coast,  is  handled  better  in  the  NRL  case  compared  to  that  using  the  NMC  analysis. 

The  SOO-mb  mass  and  wind  fields  are  very  similar  from  both  cases  at  25/1 8Z. 
However,  differences  do  exist  in  the  500-mb  vertical  motion  (Fig.  4.13)  and  moisture 
fields  (not  shown).  In  the  NRL  case  4.3,  a  significant  upward  vertical  motion  had 
developed  off  the  Carolina  coast.  These  vertical  velocities  are  missing  in  the  NMC  case 
4.4.  As  a  result,  no  precipitation  was  forecast  along  the  coastal  front  during  the  first  6  h 
of  integration  in  the  case  originating  from  the  NMC  analysis.  Case  4.3  also  produced 
stronger  vertical  motions  over  the  Mississippi  Gulf  coast  associated  with  the  secondary 
jet  circulation  and  an  area  of  positive  vorticity  advection  at  500  mb.  Both  25/12Z 
analyses,  NMC  and  NRL,  had  a  strong  500-mb  moisture  tongue  off  the  Carolina  coast. 
Unlike  that  of  the  NMC  case,  the  stronger  coastal  front  circulation  from  the  NRL  case 
had  pushed  the  moisture  tongue  inland  by  25/1 8Z. 

At  26/OOZ,  the  coastal  front  has  developed  from  Georgia  to  southern  New  England 
in  an  onshore  easterly  flow  regime.  Figure  4.14a  illustrates  that  NRL  case  still  retains  a 
stronger  thermal  gradient  associated  with  the  front  than  does  the  NMC  case  4.4  (Fig. 
4.14b).  The  inverted  surface  trough  of  low  pressure  remains  more  pronounced  in  the 
NRL  case  as  compared  to  the  NMC  case.  In  both  cases,  the  surface  low  along  the 
Alabama  coast  is  in  good  agreement  with  the  analyses  of  Bosart  and  Doyle  and  Warner 
for  position  and  intensity. 


Figure  4.1 1.  Six-hour  forecast  of  sea-level  pressure  (mb),  surface  air  temperature  (*Q, 
and  10(X)  mb  winds  valid  at  2S/18Z  for  (a)  case  4.4  from  the  NMC  analysis  and  (b)  case 
4.3  from  the  NRL  analysis.  Contours  as  in  Fig.  4.3. 
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Figure  4.12.  GALE  analysis  of  sea-level  pressure  deviation  (from  1(X)0  mb),  siuface 
temperature  (*C)  and  wind  for  25/1 8Z  by  Bosart  (1988).  Surface  temperature  contours 
every  5*C. 


Hgure  4.13.  Six-hour  forecast  of  vertical  velocity  (mb  h'*)  at  5(X)  mb  valid  at  25/1 8Z  for 
(a)  case  4.3  from  the  NRL  analysis  and  (b)  case  4.4  from  the  NMC  analysis.  Contours 
of  vertical  motion  as  in  Fig.  4.7. 
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Figure  4.14.  Twelve-hour  forecast  of  sea-level  pressure  (mb),  temperature  (*C)  and 
winds  at  10(X)  mb  valid  at  26/OOZ  for  (a)  case  4.3  from  the  NRL  analysis  and  (b)  case 
4.4  from  the  NMC  uudysis.  (c)  and  (d)  as  in  (a)  and  (b)  but  for  the  forecast  of  six-hour 
accumulated  precipitatitm  (25/1 8Z-  26/OOZ)  in  cm.  Contours  of  precipitation  at  0.25  and 
every  0.5  cm  for  0.5  and  above. 
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The  500-mb  relative  humidity  fields  are  now  very  similar;  after  12  h  of  integration, 
model  spinup  in  the  NMC  case  is  essentially  complete.  Still,  the  500-mb  absolute 
vorticity  in  the  NRL  case  over  Georgia  associated  with  the  surface  cyclone  along  the 
central  Gulf  coast  is  much  stronger  than  that  in  the  NMC  case.  In  addition,  the  upward 
vertical  motion  associated  with  the  inverted  surface  trough,  i.e.  the  southern  extent  of  the 
coastal  front,  over  the  Georgia  coast  is  much  better  developed  in  the  forecast  originating 
from  the  NRL  multivariate  analysis.  Figures  4.14c  and  d  depict  the  accumulated  6-h 
rainfall  ending  at  26/OOZ  for  the  NRL  and  NMC  cases.  The  key  difference  is  again  that 
the  NRL  case  is  able  to  produce  precipitation  due  to  mesoscale  circulations  associated 
with  the  coastal  front  along  the  Georgia  and  South  Carolina  coasts;  whereas,  the  NMC 
case  fails  to  develop  any  rainfall  in  this  region.  The  NRL  case  is  in  more  agreement  with 
the  observed  precipitation  shown  in  Fig.  4.15,  which  does  show  rainfall  associated  with 
the  coastal  front  and  with  the  low  along  the  Gulf  coast. 

The  effects  of  starting  with  stronger  secondary  circulations  in  the  entrance  region  to 
the  jet  and  in  the  coastal  front  are  illustrated  in  cross  sections  of  the  ageostrophic  flow  and 
the  relative  humidity  across  the  Carolina  coast  in  Fig.  4.16.  The  main  difference  is  the 
stronger  secondary  circulation  associated  with  the  jet  entrance  region  in  the  NRL  case.  In 
the  NRL  case  (Fig.  4.16a),  the  coastal  circulation  is  clearly  captured  by  the  jet's 
circulation,  which  also  intensifies  the  subsidence  in  the  cutoff  low  offshore  and  the 
ageostrophic  flow  across  the  coastal  front.  Higher  low-level  humidities  are  also  advected 
over  the  front  in  the  NRL  case  (Fig.  4.16b)  to  produce  the  forecast  precipitation. 

The  cyclonic  circulation  which  formed  near  Jacksonville,  Florida  continues  up  the 
coast  as  a  weak  low.  By  26/06Z,  it  iS  located  near  Wilmington,  North  Carolina.  The 
low  along  the  central  Gulf  coast  is  just  a  few  miles  southeast  of  Pensacola,  Florida.  The 
positions  of  these  circulations  from  the  NRL  case  are  in  excellent  agreement  with 
Bosart's  analysis  (not  shown).  Once  again,  the  5()0-mb  ascending  vertical  velocities 
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Figure  4.16.  Twelve-hour  forecast  of  (a)  the  vertical  motion  (mb  h*^)  and  ageostrophic 
wind  (m  s**)  and  (b)  relative  humidity  (%)  in  the  plane  of  cross  section  CS2  valid  at 
26/OOZ  for  case  4.3  from  the  NRL  analysis,  (c)  and  (d)  as  in  (a)  and  (b)  but  for  case  4.4 
from  die  NMC  analysis.  See  Fig.  4.5c  for  location  of  CS2. 
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along  the  coastal  front  and  associated  with  the  surface  wave  over  the  Florida  panhandle 
are  about  25%  greater  in  the  NRL  case  versus  the  NMC  case. 

4.4.3  Comparisons  using  26/OOZ  and  26/12Z  Cases 

The  analyses  for  26/OOZ  and  the  subsequent  forecasts  are  compared  in  cases  4.5 
and  4.6.  At  26/OOZ,  the  NRL  analysis  once  again  produced  a  much  tighter  thermal 
gradient  across  the  coastal  front  than  shown  by  the  coarser  NMC  analysis.  In  this  case 
however,  the  NRL  analysis  was  able  to  significantly  sharpen  the  first-guess  temperature 
gradient  across  the  Carolina  coast.  The  cold  air  damming  is  missing  in  the  first-guess 
field  (Fig.  4.17a)  and  NMC  analysis  (not  shown),  but  is  still  evident  in  the  NRl  analysis 
(Fig.  4.17b).  The  NRL/NCSU  analysis  scheme  was  able  to  use  the  observational  data  to 
enhance  the  background's  weaker  gradient  generated  by  our  mesoscale  model. 

At  500  mb,  the  NRL  analysis  in  case  4.5,  like  in  case  4.3,  is  successful  in 
correcting  the  background  fields.  The  phase  of  the  temperature  field  is  again  adjusted  and 
the  cutoff  low  has  been  deepened  by  48  m  and  is  within  6  m  of  the  NMC  analyzed  value. 
The  strong  vorticity  maximum  at  the  base  of  the  large-scale  trough  at  500  mb  is  of  the 
same  magnitude  in  both  cases.  A  weak  vorticity  maximum  is  present  over  Georgia  in  the 
NRL  analysis  (Fig  4.17d)  which  is  not  discernible  in  the  NMC  analysis  (Fig  4.17c). 
This  is  a  reflection  of  the  stronger  secondary  circulation  in  the  entrance  region  of  the  jet 
over  the  eastern  U.S.  as  it  interacts  with  the  coastal  low.  At  250  mb,  this  jet  is  also 
stronger  and  sharper  in  the  NRL  analysis  (not  shown).  As  the  short  wave  over  the  Great 
Lakes  moves  eastward,  it  produces  confluence  of  the  flow  into  the  jet  thereby 
strengthening  it.  The  secondary  circulation  in  the  entrance  region  of  the  jet  has  begun  to 
absorb  and  deepen  the  coastal  front's  ageostrophic  circulation  as  shown  in  the  cross 
sections  (CS2)  in  Fig.  4.18.  However,  the  coastal  front  circulation  is  still  a  distinct 
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Figure  4.17.  Sea-level  pressure  (mb)  and  1000-mb  temperature  (*C)  and  winds  at 
IfHoOZ  for  (a)  the  first  guess,  a  12-h  forecast  from  case  4.4  and  (b)  the  NRL  analysis. 
Absolute  voiticity  (in  units  of  lO"^  s'*)  with  geopotential  height  (m)  at  500  mb  at  the 
same  time  for  (c)  the  NMC  analysis  and  (d)  the  NRL  analysis. 


Figure  4.18.  The  magnitude  of  (a)  the  wind  (m  s**)  normal  to  the  plane  of  the  cross 
section  and  (b)  the  relative  humidity  (%)  in  the  cross  section  CS2  for  the  NRL  analysis  at 
26/QCIL.  The  vertical  motion  and  ageostrophic  wind  in  the  plane  of  cross  section  C^2  at 
26/(X)Z  for  (c)  the  NRL  analysis  and  (d)  for  the  NMC  analysis. 
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feature  in  the  cross  section,  particularly  for  the  NRL  analysis  shown  in  Fig.  4.18c.  The 
ageostrophic  circulation  due  to  the  coastal  front  remains  stronger  and  shaqjer  in  the  NRL 
analysis  as  compared  to  the  NMC  analysis  in  Fig.  4.18d.  Higher  humidities  are  also 
found  over  the  region  of  the  coastal  front  in  the  NRL  analysis  (Fig.  4.18b).  The 
mesoscale  circulation  derived  from  our  analysis  is  also  consistent  with  the  12-h  forecast 
produced  in  case  4.3,  giving  confidence  in  our  diagnosed  secondary  circulations.  Also, 
the  horizontal  temperature  gradient  in  the  lowest  100  mb  of  the  atmosphere  in  the  vicinity 
of  the  front  (1200  - 1600  km  from  northernmost  end  of  CS2)  is  much  greater  in  case  4.3. 

By  26/06Z,  the  forecast  from  the  NRL  case  4.5  possessed  stronger  500-mb 
vorticity  and  vertical  velocities  along  the  east  coast  than  in  the  NMC  case  4.6.  The 
resulting  6-h  precipitation  pattern  is  much  better  in  case  4.5  (Fig.  4.19a).  With  a  fine- 
resolution  model  forecast  serving  as  the  first  guess,  the  NRL  analysis  has  the  advantage 
of  possessing  the  strong  mesoscale  circulations  and  higher  initial  humidities  which  help 
to  force  the  precipitation.  Whereas,  the  forecast  in  case  4.6  (Fig,  4.19b)  originates  from 
the  coarse  NMC  analysis,  which  fails  to  capture  these  features  and  leads  to  a  substantial 
model  spin-up  which  is  readily  apparent  in  the  precipitation  pattern,  when  compared  to 
the  observed  rainfall  in  Fig.  4.19c. 

After  12  h  of  model  integration,  the  surface  temperature  gradient  in  the  coastal 
trough  is  still  markedly  stronger,  and  heavier  precipitation  is  generated  in  the  lee  of  the 
Appalachians  in  the  forecast  originating  from  the  NRL  analysis  (Fig.  4.20).  However, 
by  this  time,  the  upper-level  features  in  the  two  cases  4.5  and  4.6  are  very  similar.  A 
cross  section  (CS3  from  84*W,  43.5*N  to  62*W,  27.5*N)  across  the  coast  in  the  region 
of  the  coastal  low  shows  a  strong  vertical  circulation  in  the  coastal  low  merged  with  the 
secondary  circulation  in  the  jet  entrance  region  (Fig.  4.21a  and  b).  A  cross  section  (CS4 
from  70*W,  53*N  to  80*W,  28*N)  along  the  coastline  cutting  through  this  jet  and  the  one 
at  the  base  of  the  large-scale  trough  is  shown  in  Fig.  4.21c  and  d.  The  coastal  low 
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Figure  4.20.  Twelve-hour  forecast  of  (a)  sea-level  pressure  (mb)  and  1000-mb 
temperature  (*Q  and  winds  and  (b)  the  geopotential  height  (m),  temperature  (*C),  and 
winds  at  250  mb  valid  at  2b/12Z  for  case  4.5  from  the  NRL  analysis  at  26/12Z.  The 
forecast  6-h  accumulated  precipitation  in  cm,  valid  at  26/12Z  for  (c)  case  4.5  from  the 
NRL  analysis  and  (d)  case  4.6  from  the  NMC  analysis. 
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Figure  4.21.  Twelve-hour  forecast  for  case  4.5  of  the  (a)  wind  (m  s'l)  normal  to  cross 
section  CS3  (see  Fig.  4.20a)  and  (b)  the  votical  motion  (mb  h'^)  and  ageostrophic  wind 
(m  s'^)  in  the  plane  of  CS3,  valid  at  26/1221  (c)  and  (d)  as  in  (a)  and  (b)  but  for  the 
cross  section  CS4  (see  Fig.  4.20b). 
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circulation  is  shown  to  the  rear  of  the  strengthening  direct  secondary  circulation  in  the 
entrance  region  of  the  northern  jet.  Ascent  to  the  south,  which  is  stronger  in  the  forecast 
from  the  NRL  analysis,  is  associated  with  the  exit  region  of  the  jet  at  the  base  of  the 
large-scale  trough.  This  ascent  leads  to  the  increased  precipitation  over  western  Georgia, 
which  is  confirmed  by  the  observed  precipitation. 

The  analyses  shown  in  Fig.  4.22  at  26/1 2Z  confirm  the  large  scale  of  the 
developing  surface  low  pressure  and  upper-air  trough,  as  the  strong  short  wave  reaches 
its  base  and  the  northern  jet  is  strengthened  by  the  short  wave  from  the  Great  Lakes.  Due 
to  this  large  scale,  the  differences  between  the  NRL  and  NMC  analyses  are  less  apparent. 
The  merging  of  the  circulations  associated  with  the  entrance  region  of  the  northern  jet  and 
the  coastal  low  is  still  better  resolved  however  by  the  NRL  analysis.  Cross  sections 
along  (CS4)  and  across  (CS3)  the  coastline  in  Fig.  4.23  show  the  merging  of  the  two 
circulations.  Two  circulations  can  be  clearly  seen  in  the  cross  section  across  the  coast  in 
the  region  of  maximum  ascent  in  the  coastal  low.  As  a  result  of  the  larger  scale  of  the 
new  developing  system,  the  subsequent  forecasts  from  the  NRL  (case  4.7)  and  NMC 
(case  4.8)  analyses  reveal  less  differences  than  previously.  The  12-  and  24-h  surface 
forecasts  for  the  NRL  case  4.7  (Fig.  4.24)  show  the  development  of  a  new  low  off 
Hatteras,  which  moves  northward  and  deepens  rapidly  to  overshadow  the  earlier  coastal 
low.  In  Fig  4.25,  cross  sections  (CSS  from  70*W,  50.5’N  to  76*W,  28’N  and  CS6 
from  68.7*W,  55.5*N  to  76.7*W,  25.5’N;  see  Fig.  4.24)  along  the  coast  show  the 
secondary  circulation  in  the  exit  region  of  the  southern  jet  catching  up  with  the  circulation 
in  the  entrance  region  of  the  northern  jet,  so  that  they  cooperate  in  strengthening  the 
ascent  in  the  rapidly  developing  low.  Such  has  been  found  in  other  cases  analyzed  by 
Uccellini  and  Kocin  (1987). 


Figure  4.22.  (a)  Sea-level  pressure  (mb)  and  l(XX)-mb  tempmture  CQ  and  winds,  (b) 
the  geopotendal  height  (m),  temperature  CQ,  and  winds  at  250  mb,  (c)  absolute  vordcity 
(in  units  of  lO*^  s**)  widi  geopotendal  height  (m)  at  5(X)  mb,  and  (d)  the  verdcal  velocity 
in  mb  h*l  at  5(X)  mb  for  the  NRL  analysis  at  2^12Z. 
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Figure  4.22.  Continued. 
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Figure  4.23.  (a)  The  wind  (m  s'^)  normal  to  the  cross  section  CS3  and  (b)  the  vertical 
motion  (mb  h‘l)  and  ageostrophic  wind  (m  s'^)  in  the  plane  of  CS3  for  the  NRL  analysis 
valid  at  26/12Z.  (c)  and  (d)  as  in  (a)  and  (b)  but  for  the  cross  section  CS4. 
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Figure  4.25.  Twelve-hour  forecast  for  case  4.7  of  the  (a)  wind  (m  s‘^)  normal  to  cross 
section  CSS  (see  Fig.  4.24a)  and  (b)  the  vertical  motion  (mb  h'^)  and  ageostrophic  wind 
(m  s'l)  in  the  plane  of  the  same  cross  section  valid  at  27/OOZ.  (c)  and  (d)  the  same  as  (a) 
and  (b)  but  for  a  24-h  forecast  for  cross  section  CS6  (see  Fig.  4.24b)  valid  at  27/12Z. 
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4.5  Summary 

During  the  last  decade,  optimal  interpolation  replaced  successive  correction 
methods  as  the  dominant  objective  analysis  technique  in  operational  weather  forecasting 
systems.  However,  with  the  advent  of  Bratseth's  (1986)  analysis  scheme,  which 
converges  to  the  optimal  interpolation  solution,  the  less  expensive  approach  of  successive 
correction  has  become  a  powerful  and  attractive  alternative  analysis  method.  The 
objective  analysis  scheme  developed  for  use  with  the  NRL/NCSU  mesoscale  model 
utilizes  the  Bratseth  scheme,  in  which  the  data  weights  are  dependent  on  the  covariance 
between  observations,  are  reduced  in  regions  of  higher  data  density,  and  include 
observational  errors. 

A  computationally  inexpensive  multivariate  methodology  for  linking  the  mass  and 
momentum  fields  has  been  devised  and  is  different  from  earlier  methods  reported  in  the 
literature.  The  key  element  is  obtaining  the  gradient  of  the  geopotential  change  from  the 
change  in  u-  and  v-  through  the  geostrophic  relation  at  each  iteration  of  the  analysis 
process. 

A  qualitative  evaluation  of  the  NRL/NCSU  objective  analysis  scheme  has  been 
accomplished.  Several  tests  were  conducted  using  a  GALE  lOP  2  data  set.  NRL 
analyses  verified  well  against  Bosart's  hand  analyses.  The  multivariate  analysis  was 
successful  in  deepening  the  cutoff  low  over  the  western  Atlantic.  The  univariate  analysis 
showed  some  improvement,  but  as  expected,  the  multivariate  analysis  proved  to  be 
superior  by  using  wind  information  to  achieve  gradients  of  geopotential  through  the 
geostrophic  relationship.  The  NRL  multivariate  analysis  was  also  successful  in  bringing 
the  first-guess  temperature  pattern  into  proper  phase. 

A  very  encouraging  finding  was  that  the  second-order  fields  of  absolute  vorticity 
and  vertical  motion  were  significantly  improved  by  the  objective  analysis.  Gradients  and 
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maxima  and  minima  were  pushed  toward  those  found  in  the  NMC  analyses.  In  addition, 
the  upper-level  jet  from  the  background  field  was  strengthened  and  tightened. 

In  the  cases  where  the  forecast  originated  from  an  NRL  analysis,  the  forecast  fields 
maintained  the  mesoscale  features  better  than  those  originating  from  an  NMC  analysis. 
The  first-guess  fields  from  a  prior  model  forecast  already  contains  those  finer  scale 
features.  The  objective  analysis  was  not  only  successful  in  retaining  those  finer  scale 
features,  but  also  was  able  to  correct  their  phases  and  magnitudes.  The  mesoscale 
circulations  associated  with  the  coastal  fi'ont  on  25-26  January  1986  were  forecast  more 
accurately  when  the  NRL  multivariate  analysis  was  incorporated. 
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5.  DIABATIC  INITIALIZATION  EXPERIMENTS 
WITH  CONVECTIVE  HEATING 

The  analyses  presented  in  the  previous  chapter  were  initialized  adiabatically  for 
model  integration.  That  is,  the  vertical  mode  initialization  procedure  did  not  include 
forcing  due  to  diabatic  heating.  Here,  convective  heating  is  accounted  for  in  the  model 
initialization. 

5.1  Introduction 

The  dominant  forcing  in  the  troposphere  within  areas  of  precipitation  is  diabatic 
heating  (convective  and  radiative  heating  and  large-scale  condensation).  The  present 
study  deals  strictly  with  diabatic  heating  from  cumulus  convection.  A  convective  heating 
term  is  incorporated  in  the  balance  condition  for  the  high  frequency  normal  modes 
(Bengtsson  1981;  Errico  and  Rasch  1988;  Kitade  1983;  Wergen  1983).  With  the 
inclusion  of  heating  corresponding  to  observed  rainfall  data  in  the  diabatic  forcing  term, 
the  initialization  can  force  the  model  state  toward  the  real  atmospheric  state.  As  a  result, 
the  diabatic  initialization  corrects  the  model's  forecast  fields.  Possibly  the  primary  benefit 
of  initialization  with  diabatic  forcing  is  that  the  spinup  problem  (the  extent  that  the  initial 
amplitude  of  the  divergent  part  of  the  wind  and  the  precipitation  are  in  error)  can  be 
minimized.  That  is,  diabatic  vertical  mode  initialization  (DVMI),  also  referred  to  as 
diabatic  nonlinear  normal  mode  initialization,  produces  an  initial  model  state  which 
already  includes  convective  heating. 

Model  spinup  appears  in  all  numerical  forecast  systems  both  in  midlatitudes  and 
tropics  (Kasahara  et  al.  1988).  It  is  one  of  the  most  serious  problems  in  numerical 
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weather  prediction  and  has  been  documented  in  research  and  operational  models.  Girad 
and  Jarraud  (1982)  and  Heckley  (1985)  documented  spinup  in  the  operational  models 
used  by  the  European  Centre  for  Medium  Range  Weather  Forecasts.  Whereas,  Miyakoda 
et  al.  (1978)  and  Donner  (1988)  illustrated  different  aspects  of  the  spinup  phenomena  in 
research  models  at  the  Geophysical  Fluid  Dynamics  Laboratory  and  National  Center  for 
Atmospheric  Research,  respectively. 

Inadequate  physical  parameterizations  and  incorrect  initial  speciEcations  of 
divergence,  moisture  and  thermal  fields  are  the  primary  causes  of  spinup  (Mohanty  et  al. 
1986).  Until  recently,  most  normal  mode  initialization  methods  were  adiabatic  in  nature; 
that  is,  the  existing  latent  heating  at  the  initial  model  time  was  excluded  leaving  out  a 
significant  heat  source  which  affected  the  adjustment  between  the  mass  and  momentum 
fields.  Lejenas  (1980)  found  that  adiabatic  normal  mode  initialization,  while  providing 
improved  vertical  velocity  profiles,  does  not  alleviate  the  spinup  problem.  Data* 
assimilation  techniques  which  incorporate  a  spun-up  background  field  in  the  analysis  of 
observed  data  have  been  developed.  At  the  start  of  each  update  cycle  of  the  assimilation, 
the  background  forecast  field  is  corrected  using  the  observed  mass  and  momentum  fields. 
During  this  process,  the  divergent  part  of  the  wind  is  lost  primarily  because  the 
divergence  has  the  same  order  of  magnitude  as  the  errors  in  the  wind  field.  This 
divergent  circulation  associated  with  upward  vetical  motion  can  be  regained  if  the  analysis 
at  each  cycle  of  the  assimilation  is  initialized  with  a  diabatic  initialization  procedure. 
Turpeinen  et  al.  (1990)  concluded  that  the  inclusion  of  observed  low-level  latent  heat 
sources  with  a  consistent  and  higher-resolution  moisture  analysis  are  critical  components 
in  reducing  the  underestimation  of  initial  divergence,  and  hence  precipitation  rates. 

Throughout  the  1980s,  NWP  models  incorporating  initial  latent  heating  based  on 
rain  rates  have  been  developed.  Turpeinen  et  al.  (1990)  summarized  the  main  features  of 
some  of  these  models.  Mohanty  et  al.  (1986)  investigated  the  inclusion  of  diabatic 
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heating  in  the  vertical  mode  initialization  (VMI)  and  found  that  the  divergent  circulation 
was  rapidly  destroyed  by  the  integrated  model  heating  if  it  differed  too  much  from  the 
observed  heating.  In  the  present  study,  we  have  merged  the  model-produced  heating, 
during  the  first  few  hours  of  integration,  with  the  heating  used  in  the  initialization. 

The  major  thrust  here  is  to  incorporate  the  latent-heat  forcing  into  the  vertical  mode 
initialization  of  the  NR17NCSU  mesoscale  model  in  an  effort  to  shorten  the  spinup  time. 
The  goal  is  also  to  demonstrate  that  VMI  with  diabatic  forcing  provides  improved 
initialized  fields  over  adiabatic  VMI.  A  diabatic  heating  term,  derived  from  'spunup' 
model  convective  rainfall  and  observed  rainfall,  has  been  included  in  the  vertical  mode 
initialization  to  better  define  the  initial  conditions  for  model  integration.  Several  model 
runs  (test  cases)  are  made  to  test  the  effectiveness  of  DVMI  as  compared  to  initialization 
excluding  the  heating  term,  or  adiabatic  VMI  (A VMI).  The  procedure  for  computing  the 
diabatic  heating  for  the  initialization  is  described,  followed  by  the  experimental  design 
and  results. 

5.2  Diabatic  Heating  in  Model  Initialization 

Until  recently  the  NRL/NCSU  model  was  initialized  exclusively  using  an  adiabatic 
vertical  mode  initialization.  The  procedure  is  detailed  by  Sashegyi  and  Madala  (1990, 
1992a)  and  is  discussed  briefly  in  Chapter  3  of  this  document.  In  their  VMI  only  the  first 
three  modes,  the  external  mode  and  first  two  internal  modes  of  the  model,  are  initialized. 
Sashegyi  and  Madala  (1990)  showed  that  the  first  two  modes  require  only  two  iterations 
of  the  scheme  to  converge.  The  third  mode  converged  after  five  iterations.  The  phase 
speeds  of  these  modes  are  much  faster  than  typical  meteorological  systems,  which  naove 
at  speeds  less  than  20  m  s'^. 
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Adiabatic  VMI  has  proven  to  be  extremely  effective  in  controlling  spurious  large- 
amplitude  gravity  wave  oscillations  due  to  initial  imbalances  between  the  mass  and  wind 
fields  in  primitive  equation  models  (Puri  and  Miller  1990;  Sashegyi  and  Madala  1990). 
However,  the  use  of  A  VMI  leads  to  a  drastic  depletion  of  the  divergent  circulation. 

Wergen  (1983  and  1987)  introduced  a  procedure  to  overcome  this  problem,  in 
which  average  diabatic  heating  is  obtained  by  integrating  the  model  for  a  few  time  steps 
prior  to  the  initialization.  The  average  heating  is  then  included  in  the  nonlinear  forcing  of 
the  iterative  initialization  process.  The  fixed  diabatic  heating  is  determined  as  the  time- 
mean  heating  generated  by  model  physics  (cumulus  parameterization  and  large-scale 
heating)  during  the  short  integration.  Unfortunately,  most  forecast  models  do  not 
produce  realistic  heating  rates  during  the  first  few  hours,  leading  to  the  spinup  problem. 
Therefore,  for  present  experiments,  the  choice  was  made  to  integrate  the  model  for  9  h 
before  accumulating  the  model  latent  heating  and  convective  rainfall  from  which  our  time- 
mean  heating  and  rainfall  rates  were  derived.  The  diabatic  VMI  (DVMI)  used  in  this 
study  is  similar  to  Wergen's,  but  it  is  original  in  that  'spun-up'  model  heating  is  applied 
in  the  initialization,  rather  than  the  heating  from  the  first  few  hours  of  integration  as  is 
usually  done.  The  procedure  involves  the  following  steps: 

-  The  model  is  integrated  for  12  h  from  adiabatically  initialized  fields.  Convective 
rainfall,  from  which  the  diabatic  heating  is  derived,  is  accumulated  from  9-12  h  of  this 
initial  integration  (control  run),  allowing  ample  time  for  model  spinup. 

"  Time  average  convective  rainfall  rates  over  the  3-h  period  are  derived  for  all  grid 
points. 

-  The  average  rain  rates  are  converted  into  vertical  profiles  of  heating  using  a 
reverse  Kuo  cumulus  parameterization.  This  vertical  distribution  of  the  latent  heating  is 
based  on  the  cumulus  convective  scheme  of  the  model.  Thus,  no  simplifying  hypotheses 
are  needed.  The  advantage  of  the  present  approach  is  an  improved  consistency  obtained 
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with  the  physical  parameterization  scheme  of  the  model.  In  future  studies,  observed 
rainfall  will  be  combined  with  our  model  rainfall  rates  in  this  procedure  to  obtain  more 
objectivity. 

"  This  three-dimensional  diabatic  heating  field  is  then  added  as  a  forcing  term  in 
the  VMI  thermodynamic  energy  balance  equation  as  the  implicit  DVMI  is  performed.  As 
a  result,  diabatic  effects  are  implicitly  included  in  the  specification  of  the  heating  field. 

The  diabatic  initialization  described  above  is  primarily  for  heating  associated  with 
convective  precipitation.  The  heating  is  nearly  zero  in  the  lower  troposphere  and  varies 
significantly  from  case  to  case  in  the  upper  troposphere.  The  first  three  modes  of  the 
model  account  for  40  to  50%  of  the  convective  heating  in  the  upper  troposphere. 
Therefore,  we  assume  that  the  method  described  here  accounts  for  a  significant  part  of  the 
divergence  field. 

The  conversion  of  rainfall  rates  into  vertical  profiles  of  convective  heating  dq^ends 
on  the  stability  of  the  environment  as  determined  by  the  cumulus  parameterization.  If 
convective  instability  occurs,  a  heating  function  is  computed  using  the  reverse  Kuo 
scheme.  This  heating  function  is  then  used  to  distribute  the  diabatic  heating,  derived 
from  the  average  rainfall  rates,  in  the  vertical  direction.  If  the  environment  is 
convectively  stable,  the  diabatic  heating  associated  with  the  two-dimensional  rainfall 
distribution  will  be  taken  as  zero.  Therefore,  non-convective  precipitation  does  not 
contribute  to  the  diabatic  heating  in  the  present  initialization  scheme. 

The  persistence  of  balance  achieved  by  DVMI  during  the  early  stages  of  the  model 
integration  depends  heavily  on  the  compatibility  between  the  specified  heating  during 
initialization  and  the  model-produced  heating  during  integration  (Puri  and  Miller  1990). 
Puri  (1987)  suggested  that  one  possible  way  of  retaining  this  compatibility  in  the  heating 
rates  would  be  to  adjust  the  moisture  field  until  the  model-produced  heating  rates 
determined  by  the  convective  parameterization  in  the  model  are  similar  to  those  used 
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during  DVMI.  Here,  two  procedures  are  used  to  retain  this  desired  compatibility  in 
heating  rates.  First,  the  model  rainfall  and  heating  rates  produced  during  the  first  3  h  of 
integration  are  merged  with  the  corresponding  rates  from  the  initialization.  A  nonlinear 
weighting  factor  (a)  is  employed  where  a=l  at  t=0  h  and  a=0  at  t=3  h.  The  relation  is 
given  by: 

Rni  =  aRi  +  (l-a)Rc  5.1 

where  Rj  is  the  convective  rainfall  rate  used  in  the  diabatic  VMI,  Rp  is  the  model- 
produced  convective  rainfall  rate,  R^  is  the  merged  rain  rate,  and  a  is  given  by: 

[l +sin(7Ct/3-7t/2)]  52 

2 

The  above  relations  allow  the  influence  of  the  initialized  heating  and  rain  rates  to  approach 
zero  in  a  smoother  manner  than  a  simple  linear  relationship  would  allow.  Merging  is 
applicable  only  for  DVMI  cases  since  AVMI,  by  definition,  does  not  include  a  diabatic 
heating  term  in  its  balance  equation. 

The  second  method  involves  adjusting  the  humidity  fields  prior  to  initialization. 
The  imponance  of  moisture  adjustment  in  diabatic  initialization  was  emphasized  by 
Wolcott  and  Warner  (1981).  They  argued  that  if  the  environment  is  not  humid  enough, 
the  upward  motion  associated  with  the  initialized  divergence  field  will  not  be  sustained  by 
latent  heat  release.  In  the  present  study,  specific  humidity  is  enhanced  only  at  grid  points 
where  convective  rainfall  is  occurring  and  the  relative  humidity  is  less  than  95%.  At 
those  points,  the  specific  humidity  is  increased  to  95%  of  the  saturation  specific  humiditj’ 
at  the  three  lowest  model  levels.  However,  our  findings  indicated  that  the  changes  in  the 
humidity  fields  were  small,  and,  as  a  result,  the  adjustment  had  minimal  affect  on  the 
convective  rainfall.  Since  the  modeled  large-scale  condensation  occurs  only  upon 
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supersaturation,  not  at  95%  relative  humidity,  this  moisture  adjustmeru  would  nave  little 
affect  on  non-convective  precipitation. 

5.3  Experimental  Design 

Initially,  three  test  cases  (model  runs)  were  performed  to  determine  improvements 
(if  any)  in  the  model  initial  state  and  subsequent  forecasts  using  DVMI  versus  AVMI. 
The  first  is  a  control  run  (case  5.1)  in  which  the  model  is  initialized  at  1200  UTC  26 
January  1986  (26/12Z)  using  the  NMC  analyzed  data.  These  data  of  2.5*  resolution  are 
first  interpolated  to  the  model  grid  and  then  initialized  with  the  AVMI  for  the  first  three 
vertical  normal  modes  of  the  model.  The  model  is  then  integrated  for  24  h.  The  12-h 
forecast  from  the  control  run  serves  as  the  initial  fields  for  cases  5.2  and  5.3.  In  case 
5.2,  the  model  is  initialized  using  DVMI  on  the  forecast  fields  from  case  5.1  valid  at 
27/OOZ.  Model-generated  rainfall  (26/21Z  to  27/OOZ  of  the  control)  is  used  to  compute 
heating  rates  at  each  vertical  level  of  the  model  for  the  DVMI.  The  model  is  then 
integrated  for  12  h,  with  heatA’ain  merging.  The  initialization  for  the  next  case  (case  5.3) 
is  simply  an  AVMI  on  the  control  run's  12-h  forecast. 

For  this  study,  the  control  run  is  considered  the  "truth"  for  comparison  against  the 
other  two  cases.  From  26/1 2Z  to  27/OOZ  of  the  control  run,  the  model  should  have 
sufficient  time  to  spin  up,  or  in  other  words,  develop  the  divergent  part  of  the  wind.  By 
using  model-generated  rainfall  rates  in  the  DVMI  (case  5.2),  one  would  expect  a 
negligible  spinup  during  the  subsequent  forecast.  Also,  one  might  expect  little  change 
due  to  the  initialization  using  model  rain  since  the  model  heating  function  is  consistent 
with  other  model  parameters.  By  comparing  the  initialized  and  subsequent  forecast  fields 
from  case  5.2  with  the  corresponding  fields  from  the  control,  one  can  determine  if  the 
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model-generated  rainfall,  which  is  a  two-dimensional  (2-D)  distribution,  can  re-create  the 
3-D  distribution  of  vertical  heating  rates  accurately. 

5.4  Results 

The  following  results  were  obtained  by  making  qualitative  comparisons  of 
meteorological  fields  firom  the  different  model  initializations  and  forecasts  of  the  three  test 
cases.  The  differences  between  cases  5.2  and  5.3  of  their  respective  initialized  fields  of 
1000  mb,  500  mb,  and  250  mb  geopotential  height,  temperature,  and  winds  are 
extremely  small.  As  expected,  the  first-order  analysis  fields  (12-h  forecast  from  the 
control  run)  changed  little  with  each  initialization  in  these  test  cases  and,  as  a  result,  the 
forecasts  (valid  at  27/1 2Z)  of  these  variables  from  these  two  runs  are  almost  identical  and 
agree  well  with  the  control  run  forecast. 

5.4. 1  Comparison  of  Initialized  Vertical  Velocity 

Since  vertical  velocity  or  omega  (dp/dt)  is  directly  related  to  diabatic  heating,  we 
can  expect  noticeable  differences  between  the  initial  states  of  omega  from  cases  5.2  and 
5.3.  Furthermore,  the  effectiveness  of  diabatic  heating  in  inducing  divergent  circulations 
can  be  seen  by  comparing  omega  fields  from  DVMI  and  those  resulting  from  AVMI.  As 
expected,  in  geographic  areas  absent  of  precipitation,  vertical  velocity  patterns  in  cases 
5.2  (DVMI)  and  5.3  (AVMI)  are  approximately  the  same  and  very  close  to  the  "truth" 
(Fig.  5.1).  Off  the  east  coast,  an  area  of  strong  ascent  stretches  along  the  eastern 
seaboard  and  is  associated  with  the  old  frontal  boundary  and  the  developing  cyclone  off 
the  North  Carolina  coast.  In  this  region  where  the  control  run  ("truth")  predicts  heavy 
rain  at  initialization  time  (27/OOZ),  the  initialized  state  for  case  5.2  is  very  close  to  the 


a 


55.0 


LONG 


Figure  5.1.  5(X)-mb  vertical  motion  fields  valid  at  27/OOZ:  (a)  from  12-h  fr>recast  of  the 
control  (case  5.1),  (b)  after  DVMI  (case  5.2),  and  (c)  after  AVMI  (case  5.3).  Vertical 
motion  of  magnitude  of  -20  mb  h'l  or  more  are  contoured  every  2.5  mb  h'^. 
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"truth."  The  pattern  and  positions  of  maxima  and  minima  in  the  omega  field  show  little 
change  from  the  control  run,  suggesting  the  3-D  strucmre  of  model  convective  heating 
can  be  re-created  from  a  2-D  rain  distribution.  A  further  comparison  between  case  5.2 
and  case  5.3  reveals  the  deficiency  of  AVMI.  In  the  vicinity  of  the  developing  cyclone 
off  Cape  Hatteras,  the  magnimde  of  omega  has  changed  substantially  using  AVMI;  the 
maximum  upward  velocity  is  only  about  60%  of  that  resulting  from  DVMI.  Since  a 
significant  amount  of  the  divergence  (and  hence  vertical  motions)  associated  with  heavy 
precipitation  events  may  be  forced  by  latent  heat  release,  and  since  the  nonlinear  step  in 
the  AVMI  does  not  include  this  diabatic  forcing,  the  divergence  due  to  this  forcing  will  be 
removed  as  if  it  were  associated  with  unwanted  gravity  waves.  As  a  result,  the  vertical 
velocities  are  underestimated. 

5,4.2  Forecast  Precipitation  Comparisons 

The  patterns  of  12-h  (27/OOZ  -  27/1 2Z)  convective  rainfall  in  cases  5.2  and  5.3  are 
very  similar  to  each  other  and  to  the  control  in  areal  coverage  and  accumulation.  The  12- 
h  convective  rainfall  forecast  from  the  adiabatically  initialized  fields  is  just  as  good  as  the 
forecast  originating  from  the  diabatically  initialized  fields.  Since  the  physics  and 
dynamics  of  the  nonlinear  normal  mode  procedure  are  internally  consistent  with  the 
numerical  model  (which  provided  the  analysis  fields),  the  initialization  procedure, 
whether  with  or  without  diabatic  heating,  changes  the  model  forecast  fields  only  a  small 
amount.  As  a  result,  the  subsequent  forecasts  arc  essentially  the  same.  Therefore, 
additional  model  experiments  were  performed  in  which  GALE  3 -hourly  surface  and 
upper-air  data  were  analyzed  using  a  multivariate,  successive  correction  objective  analysis 
scheme  (Sashegyi  et  al.  1992).  Model  forecast  fields  served  as  the  first  guess  for  the 
analyses.  Adiabatic  (AVMI)  and  diabatic  (DVMI)  initializations  were  performed  on  the 
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analysis  fields.  In  the  latter  case,  observed  rain  data  was  combined  with  the  model- 
produced  rainfall  for  computation  of  the  3-dimensional  heating  function  in  the  model 
initialization  procedure. 

Figure  5.2  shows  6-h  observed  rainfall  and  accumulated  forecast  rainfall  valid 
26/(X)Z  -  26/06Z.  The  forecast  in  Fig.  5.2a  was  initialized  with  DVMI  at  26/OOZ.  The 
rainfall  pattern  is  similar  to  the  observed;  in  particular,  over  Alabama  and  along  the  central 
Gulf  coast  region,  the  model  was  able  to  predict  the  precipitation  during  this  short-range 
forecast.  However,  the  forecast  originating  from  AVMI  (Fig.  5.2b)  fails  to  capture  any 
rainfall  over  the  Gulf  coast  associated  with  the  developing  short  wave.  The  absence  of 
diabatic  heating  from  the  VMI  has  caused  a  spinup  problem  in  this  case.  RMS 
differences  between  these  forecasts  and  observed  rainfall  were  computed.  The  differences 
were  computed  only  at  gridpoints  in  which  at  least  0.01  cm  of  either  forecast  or  observed 
rainfall  accumulated  during  the  6-h  period.  When  initializing  with  AVMI  and  DVMI,  the 
rms  differences  were  0.9  cm  and  0.7  cm,  respectively.  The  larger  rms  difference  in  the 
AVMI  case  is  due  to  model  spinup.  The  disparity  between  these  two  cases  would  be 
even  greater,  if  only  the  accumulated  rainfall  from  the  first  3  h  of  model  integration  was 
compared.  A  complete  description  of  these  experiments  and  discussion  of  results  within 
the  context  of  four-dimensional  data  assimilation  is  given  by  Harms  et  al.  (1992c) 

5.5  Summary 

In  this  study,  the  viewpoint  taken  was  that  the  "best"  initialized  states  were  those 
which  were  least  altered  by  initialization.  From  comparisons  of  vertical  velocity  fields, 
we  found  that  DVMI  with  model  rain  produced  an  initial  state  very  similar  to  the  "truth" 
and,  thus,  successfully  re-created  the  3-D  disttibution  of  convective  heating  from  the  2-D 
model  rain  distribution.  The  AVMI  (case  5.3)  changed  the  "true"  analysis  significantly 
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Figure  5.2.  Forecasts  of  accumulated  total  precipitation  in  cm  (26A)0Z  -  26A)6Z) 
originating  firon  (a)  DV^  and  (b)  A  VML  (c)  Observed  precipitation  in  cm  for  the  same 
period.  Contours  at  0.25  cm  and  every  0.5  cm  fev  0.5  and  above. 
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by  reducing  the  maximum  vertical  velocities  to  about  60%  of  that  resulting  from  DVMI. 
With  the  diabatic  heating  term  which  forces,  to  a  large  degree,  convective  precipitation, 
the  DVMI  (unlike  the  AVMI)  retains  the  divergence  (or  vertical  motion)  associated  with 
the  convection. 

When  DVMI  is  performed  on  the  analysis  fields,  model  spinup  time  is  reduced. 
The  pre-forecast  integration  of  the  model  ensures  that  the  analyses  are  spunup  and  that  the 
circulations  in  these  analyses  are  consistent  with  the  model  resolution  and  external 
parameters  such  as  topography.  The  initializations  and  forecasts  were  performed  in  this 
study  for  the  midlatitudes,  which  dictates  a  relatively  short  spinup.  In  the  data-sparse 
tropics,  spinup  can  be  expected  to  be  much  greater.  In  the  tropical  case,  DVMI  with 
merging  will  be  critical  in  reducing  spinup  time. 
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6.  IMPACT  OF  ASSIMILATING  3-HOURLY  GALE  DATA 

The  discussion  to  this  point  has  dealt  with  the  major  components  of  data 
assimilation:  objective  analysis,  initialization,  and  model  forecast.  In  the  following 
discussion,  these  components  are  linked  together  to  form  a  four-dimensional  data- 
assimilation  system. 

6.1  Introduction 

A  limited-area  intermittent  data-assimilation  system  is  used  to  assimilate  3-hourly 
data  collected  during  the  Genesis  of  Atlantic  Lows  Experiment  (GALE),  which  was 
conducted  over  the  southeastern  U.S.  during  the  winter  of  1986.  The  system,  developed 
jointly  by  NRL  and  NCSU,  is  used  to  generate  initial  conditions  for  the  NRL/NCSU 
mesoscale  model.  The  NRIVNCSU  data-assimilation  system  consists  of  a  multivariate, 
successive  correction  objective  analysis,  a  diabatic  vertical  mode  initialization  (VMI),  and 
a  mesoscale  numerical  forecast.  The  development  of  four-dimensional  data-assimilation 
systems  for  limited  areas  is  being  driven  by  a  need  to  optimally  utilize  data  from  the  new 
observing  systems  of  the  1990s  (e.g..  Next  Generation  Radar  [NEXRAD],  Special 
Sensor  Microwave/  Imager  [SSM/I],  and  Wind  Profilers)  in  mesoscale  models.  Until 
these  new  observing  systems  become  operational,  other  data  sources  must  be  used  to  test 
assimilation  techniques. 

During  intensive  observing  periods  (lOPs)  of  GALE,  3-hourly  upper-air 
soundings  were  taken  at  selected  National  Weather  Service  (NWS)  sites  over  the  eastern 
U.S.,  supplemented  by  more  frequent  and  higher  density  soundings  taken  at  sites  along 
the  east  coast  and  by  dropsondes  deployed  offshore  by  aircraft  (Dirks  et  al.  1988).  These 
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3-hourly  asynoptic  upper-air  observations  from  the  GALE  experiment  are  ideal  for  testing 
the  impact  of  assimilating  high-frequency  observations  in  a  mesoscale  model. 

The  influence  of  GALE  data  on  the  regional  analyses  and  assimilation  system  at  the 
National  Meteorological  Center  (NMC)  has  been  investigated  by  Rodgers  et  al.  (1990) 
and  DiMego  et  al.  (1989).  Using  the  higher  density  GALE  soundings  along  the  east 
coast  at  the  synoptic  hours  of  0000  and  1200  UTC  only,  Rodgers  et  al.  (1990)  found 
only  a  minimal  impact  on  the  NMC  operational  regional  analyses  and  the  subsequent 
forecasts  which  were  rerun  for  lOPs  1  and  2.  Some  noise  and  bias  problems  were  found 
however  in  the  geopotential  in  the  mid-troposphere  in  the  supplemental  soundings  along 
the  east  coast,  which  may  have  contributed  to  the  negative  result.  Also,  at  the  time  of 
their  study,  a  6-h  forecast  from  the  global  spectral  model  was  being  used  as  the  first 
guess  for  the  operational  regional  analyses  (Hoke  et  al.  1989)  and  no  feedback  between 
the  regional  model  and  the  analysis  could  take  place.  To  test  an  intermittent  assimilation 
for  use  with  future  high-frequency  wind  profiler  data,  DiMego  et  al.  (1989)  used  a  3-h 
update  cycle,  with  the  regional  model  providing  the  first  guess,  to  assimilate  the  3-hourly 
GALE  asynoptic  data  from  a  single  12-h  period  during  lOP  1.  For  their  single  case,  an 
improvement  in  the  analysis  at  the  end  of  the  12-h  period  of  assimilation  was 
demonstrated.  A  stronger  jet  with  stronger  horizontal  and  vertical  shears  in  the  final 
analysis  led  to  a  better  definition  of  the  indirect  circulation  in  the  exit  region  of  the  jet  in 
the  initialized  field.  They  also  asserted  that  the  low-level  cyclogenesis  was  improved  in 
the  first  12  h  of  the  subsequent  forecast.  Even  assimilating  the  winds  alone  improved  the 
final  analysis  when  compared  to  the  operational  analysis  which  did  not  use  the  GALE 
asynoptic  data. 

In  the  present  paper,  a  finer-resolution  mesoscale  model  is  used  in  the  assimilation 
experiments  with  the  GALE  asynoptic  data,  and  the  model  and  analysis  domains  are 
centered  on  the  GALE  region.  The  horizontal  resolution  of  the  fine  grid  in  the 
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NRL/NCSU  model  is  about  50  km  at  45*N,  whereas  the  fine  mesh  of  NMCs  nested  grid 
model  used  by  DiMego  et  al.  (1990)  had  at  that  time  a  resolution  of  about  85  km  at  the 
same  latitude  (Hoke  et  al.  1989).  In  the  present  case,  the  data  assimilation  is  extended  for 
a  full  2  1/2-day  period  for  the  second  lOP  of  GALE,  during  which  time  a  substantial 
number  of  asynoptic  observations  were  available.  During  this  time,  several  upper 
troposphere  jets  propagated  across  the  region  of  the  GALE  data,  and  a  cyclone  developed 
in  the  region  along  the  east  coast.  Since  the  ultimate  test  of  any  data-assimilation 
system's  usefulness  is  the  quality  of  the  forecast  produced  from  the  assimilated  states, 
focus  is  placed  heavily  on  comparisons  of  subsequent  forecasts  as  well  as  on  the  analyses 
from  which  they  were  initialized. 

As  noted  by  Kanamitsu  (1989),  a  number  of  advantages  are  realized  when  using  a 
four-dimensional  data  assimilation  system  (such  as  the  one  described  in  this  paper)  for 
specifying  the  initial  state  for  a  numerical  forecast:  1)  Using  a  comprehensive 
assimilation  model  with  full  physics  as  the  assimilation  nxxlel  provides  the  most  accurate 
first  guess  for  the  following  objective  analysis.  This  allows  optimum  use  of  model 
feedback  in  regions  where  data  are  sparse  in  time  and  space.  Plus,  it  ensures  that  the 
initial  state  specification  possesses  a  resolution  compatible  with  the  high-resolution 
forecast  model.  2)  Using  the  model  forecast  as  the  first  guess  for  the  next  analysis 
enables  information  from  earlier  observations  to  be  carried  forward  to  provide  an 
independent  source  of  information  to  be  added  to  the  newly-acquired  data.  That  is,  the 
background  forecast  feeds  information  forward  to  the  next  forecast.  3)  Using  a 
multivariate  objective  analysis  ensures  a  degree  of  balance  between  the  mass  and 
momentum  fields,  which  helps  minimize  changes  due  to  the  initialization.  4)  A 
multivariate  analysis  scheme  also  enables  one  model  variable  to  be  updated  from 
observations  of  another  variable. 
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The  primary  purpose  here  is  three-fold:  1)  to  describe  our  data  assimilation  system, 
2)  to  evaluate  the  effectiveness  of  the  system,  and  3)  to  assess  what  impact  the  high- 
density  3-hourly  GALE  data  have  on  the  NRL  analyses  and  forecasts.  More  specifically, 
the  analysis  performance  is  addressed  both  by  data-fitting  statistics  and  by  comparison 
with  operational  analyses.  Also,  precipitation  forecasts  are  compared  using  predictions 
originating  from  NMC  operational  analyses,  from  assimilated  fields  excluding  the  3- 
hourly  data,  and  from  assimilated  analyses  including  the  GALE  asynoptic  data. 

6.2  The  Data-Assimilation  System 

The  NRL/NCSU  four-dimensional  data-assimilation  system  uses  the  intermittent 
data-assimilation  method  (or  analysis/fotecast  cycle)  in  which  a  nnodel  forecast  serves  as 
the  background  field,  or  first  guess,  for  the  analysis  of  new  data.  The  assimilation 
system  analyzes  asynoptic  GALE  surface  and  upper-air  data  every  3, 6,  or  12  h  using  the 
NRL/NCSU  multivariate,  successive  correction  scheme.  The  analyzed  mass  and 
momentum  fields  are  balanced  by  an  initialization  procedure.  A  short-range  forecast  is 
then  made  to  serve  as  the  first  guess  for  the  next  analysis.  The  NRL/NCSU  mesoscale 
model  is  used  as  the  assimilation  model  to  generate  the  background  fields  for  the  analysis 
at  each  cycle  of  the  assimilation.  The  same  model  is  used  to  produce  subsequent 
forecasts  from  the  assimilated  states,  thus  ensuring  compatibility  of  model  resolution  and 
physical  processes  between  assimilated  and  forecast  states.  In  addition.  Bell  (1986) 
documented  that  model  spinup  time  can  be  substantially  reduced  if  the  forecast  model 
being  initialized  is  also  the  assimilation  vehicle.  He  recognized  that  the  analyzed  fields 
from  the  assimilation  must  contain  detail  compatible  with  the  resolution  of  the  forecast 
model  grid;  otherwise,  the  model  will  require  time  to  generate  the  small  scale  circulations. 
Furthermore,  with  a  fine  grid  resolution  of  approximately  50  km,  the  NRL/NCSU 
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mesoscale  model  is  designed  to  capture  and  retain  the  mesoscale  information  from  the 
GALE  data  network. 

The  NRL/NCSU  data-assimilation  system  is  also  fully  diabatic  incorporating  3- 
hourly  observed  and  model-generated  precipitation  rates  in  its  vertical  mode  initialization. 
The  initialization  scheme  is  used  to  correct  the  divergent  part  of  the  wind  produced  by  the 
analysis  and  to  provide  initial  conditions  for  model  integration  at  the  next  assimilation 
cycle. 

At  the  start  of  each  update  cycle  of  the  assimilation,  the  analysis  procedure  corrects 
the  background  forecast  field  using  the  observed  mass  and  momentum  fields.  During 
this  process,  the  divergent  part  of  the  wind  is  lost  primarily  because  the  divergence  has 
the  same  order  of  magnitude  as  does  the  errors  in  the  wind  field.  This  divergent 
circulation  )n  be  regained  If  the  analysis  at  each  cycle  of  the  assimilation  is  initialized 
with  a  diabatic  initialization  procedure.  In  the  present  assimilation  system,  the 
interpolated  analyses  are  initialized  separately  on  both  coarse  and  fine  grids  for  the  first 
three  vertical  modes  of  the  numerical  model  using  the  vertical  mode  scheme  of  Bourke 
and  McGregor  (1983),  as  applied  to  the  NRL/NCSU  model  by  Sashegyi  and  Madala 
(1992a)  and  Harms  et  al.  (1992a).  This  scheme  generates  a  balanced  divergent  pan  of 
the  wind  field,  while  producing  only  small  changes  to  the  geopotential  and  vorticity. 

Harms  et  al.  (1992a)  discussed  the  importance  of  including  diabatic  forcing  in  the 
model  initialization.  With  the  inclusion  of  heating  corresponding  to  observed  rainfall  in 
the  diabatic  forcing  term  of  the  thermodynamic  balance  equation,  the  initialization  can 
force  the  model  state  toward  the  real  atmospheric  state  and  substantially  reduce  model 
spinup.  In  our  diabatic  vertical  mode  initialization  procedure,  diabatic  heating  rates  are 
computed  from  a  merged  field  of  observed  and  model-produced  rainfall. 

Model  convective  rainfall  is  accumulated  during  the  integration  of  the  previous 
assimilation  cycle  from  3  h  before  the  current  initialization  time  to  3  h  past  that  initial 
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time.  In  other  words,  the  model  convective  precipitation  is  accumulated  for  6  h  centered 
around  the  beginning  of  the  current  assimilation  cycle.  As  a  result,  the  model  integration 
for  each  assimilation  cycle  is  3  h  longer  than  the  cycle  duration  (period  between  updates). 
For  cycle  durations  of  3,  6,  and  12  h,  the  model  is  integrated  for  0,  3,  and  9  h, 
respectively,  before  the  rainfall  accumulation  begins.  Therefore,  in  the  latter  two 
assimilation  cases,  the  model  is  spunup  prior  to  collecting  the  model  convective  rain. 
During  the  vertical  mode  initialization,  this  model-produced  rainfall  is  combined  with 
analyzed  6-h  observed  rainfall  (satellite-derived  and  GALE  3-hourly  raingauge  data, 
described  in  the  next  section)  valid  for  the  same  period,  a  time  frame  extending  to  3  h 
beyond  the  current  analysis  time.  This  procedure  requires  the  assimilation  to  occur  at 
least  3  h  after  real  time.  However,  in  tliis  sense,  diabatic  heating  rates  used  in  the 
initialization  are  derived  from  rainfall  accumulations  which  are  time-centered  at  the 
appropriate  time  (valid  at  initialization)  and,  therefore,  are  temporally  representative  rates. 
The  observed  values  are  used  where  available  over  the  model  grid;  otherwise,  the  model 
values  are  used. 

This  6-h  combined  rainfall  field  is  converted  into  time-average  rainfall  rates  for  all 
grid  points.  The  average  rain  rates  are  converted  into  vertical  profiles  of  heating  using  a 
reverse  Kuo  cumulus  parameterization.  This  vertical  distribution  of  the  latent  heating  is 
based  on  the  cumulus  convective  scheme  of  the  model.  Thus,  no  simplifying  hypotheses 
are  needed.  The  advantage  of  the  present  approach  is  an  improved  consistency  obtained 
with  the  physical  parameterization  scheme  of  the  model.  This  three-dimensional  diabatic 
heating  field  is  then  added  as  a  fixed  forcing  term  in  the  thermodynantic  energy  balance 
equation  as  each  iteration  of  the  vertical  mode  initialization  is  performed. 

During  the  first  3  h  of  model  integration  of  each  assimilation  cycle,  the  model- 
produced  convective  rainfall  and  heating  rates  are  merged  with  the  corresponding  rates 
used  in  the  previous  initialization,  thus  ensuring  the  desired  compatibility  with  the  heating 
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rates  at  initialization.  As  mentioned  earlier,  the  heating  rates  used  in  the  initialization  are 
derived  from  a  6-h  accumulation  of  rainfall  centered  on  the  initialization  time.  Therefore, 
this  6-h  amount  includes  rainfall  accumulated  for  3  h  beyond  the  valid  time  of  the 
initialization,  that  is,  the  same  3  h  of  subsequent  model  integration  in  which  the  merging 
of  heat/rain  rates  occurs.  Thus,  these  heating  rates  used  in  the  initialization  and  first  3  h 
of  subsequent  "assimilation"  model  integration  are  representative  of  "observed"  heating 
rates  for  this  period.  Thus,  the  merging  process  forces  the  subsequent  model  heating 
toward  the  "observed"  heating. 

6.3  Experimental  Design 

The  system  was  tested  using  data  collected  over  a  2  1/2-day  period  from  1200 
UTC  25  Januaiy  to  0000  UTC  28  January  1986  (25/1 2Z  to  28/OOZ).  This  data  included 
a  set  of  upper  air  soundings,  which  were  interpolated  to  10-mb  levels  from  the  surface  up 
to  100  mb,  and  surface  observations  for  lOP  2  of  GALE.  These  datasets  were  obtained 
from  the  GALE  Data  Center  (GDC)  at  Drexel  University  (Mercer  1987).  The  data  domain 
included  the  eastern  U.S.  and  western  North  Atlantic  from  1 15*W  to  45' W  in  longitude 
and  from  lO'N  to  60’N  in  latitude.  A  higher  density  of  upper-air  and  surface  data  is 
found  in  the  inner  GALE  region  covering  the  Carolinas  and  southern  Virginia  from  the 
Gulf  Stream  to  the  Appalachians.  The  sounding  network  for  GALE  consisted  of  39 
National  Weather  Service  (NWS)  sites  which  made  launches  at  3-h  intervals  when 
weather  events  dictated.  These  observing  sites  were  supplemented  by  seven  Cross-chain 
LORAN-C  Atmospheric  Sounding  System  (CLASS)  sites,  which  made  launches  as  often 
as  every  90  minutes,  in  North  Carolina  and  South  Carolina.  Plus,  numerous 
dropwindsondes  were  deployed  over  the  Atlantic  Ocean  and  the  Gulf  of  Mexico  during 
lOP  2,  as  well  as  observing  systems  aboard  two  ships  in  Atlantic  coastal  waters  (Mercer 
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and  Kreitzberg  1986).  The  GALE  datasets  included  upper-air  observations  (primarily 
synoptic)  from  several  other  NWS  sites,  the  Canadian  Weather  Service  sites,  and  a  few 
sites  in  Mexico.  The  GALE  observations  were  grouped  into  data  windows  of  3  h 
centered  on  the  analysis  times  of  0000, 0300,  0600,  0900,  1200,  1500,  1800,  and  2100 
UTC. 

As  mentioned  in  Chapter  4,  in  regions  of  no  data  over  the  ocean,  bogus  data, 
derived  from  the  difference  between  the  first  guess  forecast  and  the  NMC  2.5* 
hemispheric  analysis,  are  used  to  merge  the  analysis  smoothly  with  the  NMC  analysis 
over  the  data-sparse  ocean  regions  at  synoptic  hours.  At  off-synoptic  times,  no  bogus 
soundings  are  added  and  the  analysis  returns  to  the  background  in  regions  of  no  data.  As 
an  example.  Fig.  6.1  shows  the  locations  of  surface  and  upper-air  data  at  26/()0Z  and 
26/03Z  after  the  quality  control.  Very  few  upper-air  real  observations  (shown  by  the 
crosses)  are  found  over  the  oceans,  where  NMC  "bogus"  data  (shown  by  stars)  have 
been  added  at  the  synoptic  hours  (Figs.  6.1a  and  c).  At  the  asynoptic  hours  (Figs.  6.1b 
and  d),  the  number  and  aerial  coverage  of  the  data  is  much  less  than  at  the  synoptic  times 
as  noted  above.  Table  6.1  shows  the  average  number  of  surface  and  upper-air 
observations,  before  and  after  quality  control,  available  at  each  analysis  time,  with  the 
number  of  bogus  data  shown  in  parentheses.  As  was  also  clear  in  Fig.  6.1,  there  is 
considerable  variation.  At  the  asynoptic  times,  only  a  limited  number  of  sounding  data  is 
available  with  these  being  solely  from  the  GALE  network.  The  actual  number  of  data 
analyzed  is  less  than  the  number  of  real  GALE  observations,  particularly  in  the  case  of 
surface  data.  In  the  quality  control  procedures,  gross  error  checks  and  buddy  checks 
eliminate  a  few  data.  Data  in  close  proximity  of  each  other  or  observations  from  the  same 
location  within  a  given  3-h  data  window  are  averaged  to  create  super  observations. 

Three-hour  observed  rainfall  amounts  from  lOP  2  of  GALE  are  used  in  the  diabatic 
vertical  mode  initialization.  The  raingauge  data  obtained  from  the  GALE  Data  Center  at 
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TABLE  6.1.  Average  number  of  observations,  before  and  after  quality  control  and 
creation  of  super  observations,  in  each  analysis  time  window.  Average  number  of  bogus 
observations  is  in  parentheses 


Analysis  Time 

Upper-air 

Saifagg 

0000  UTC 

101 

89 

(93) 

337 

170 

(64) 

0300  UTC 

49 

42 

226 

85 

0600  UTC 

51 

44 

277 

120 

0900  UTC 

48 

42 

232 

87 

1200  UTC 

113 

103 

(76) 

376 

198 

(42) 

1500  UTC 

47 

40 

232 

87 

1800  UTC 

49 

42 

308 

138 

2100  UTC 

41 

35 

225 

81 
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Drexel  University  were  bilinearly  interpolated  to  the  fine  grid  of  the  NRL/NCSU  model. 
Values  for  the  coarse  grid  were  then  obtained  by  computing  a  9-point  average  using  fine 
grid  values.  These  raingauge  data  were  supplemented,  or  merged,  with  the  satellite- 
derived  rainfall  analyses  provided  by  the  University  of  Wisconsin  (Martin  et  al.  1988). 
These  satellite-derived  rainfall  values  were  available  as  6-h  accumulations  on  a  one  degree 
latitude  by  one  degree  longitude  grid  which  extended  from  approximately  25*N  to  46'N 
and  60*W  to  85‘W.  The  6-h  derived  amounts  were  temporally  interpolated  to  obtain  3-h 
accumulations,  which  were  then  bilinearly  interpolated  to  the  model  grids.  These 
interpolated  fields  of  satellite-derived  rainfall  were  merged  with  the  GALE  raingauge  data 
to  produce  our  "observed"  3-h  rainfall  data.  When  and  where  available,  the  raingauge 
values  were  kept;  otherwise,  the  merged  values  were  represented  by  the  satellite-derived 
rainfall. 

To  test  the  effectiveness  of  the  NRITNCSU  data-assimilation  system  and  to  assess 
the  impact,  if  any,  of  the  3-hourly  GALE  data  on  the  assimilation  and  subsequent 
forecasts  initialized  from  our  assimilated  states,  five  experimental  cases  were  performed, 
as  outlined  in  Table  6.2.  In  experiment  1, 12-hourly  National  Meteorological  Center  2.5’ 
hemispheric  analyses  (NMC  analyses)  for  the  period  25/12Z  to  28/OOZ  are  used.  No  re- 
analysis  was  performed  to  incorporate  the  high-density  GALE  data;  therefore,  these 
NMC  analyses  represent  our  "no  assimilation"  case  and  serve  as  the  control  experiment 
for  comparison  with  analyses  from  the  assimilation  experiments.  Differences  in  the 
assimilated  analyses  from  the  control  analyses  are  attributable  to  both  the  assimilation  of 
the  3-hourly  data  and  the  use  of  the  high  horizontal  resolution  NRIVNCSU  mesoscale 
model  in  the  assimilation  process. 

The  25/()0Z  NMC  analysis  was  used  to  cold  start  the  remaining  experiments.  The 
first  analysis  of  GALE  data  occurs  at  25/12Z  in  these  assimilation  runs.  The  25/()0Z 
NMC  analysis  was  initialized  and  then  the  model  was  integrated  for  12  h.  The  12-h 
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TABLE  6.2.  Data  Assimilation  Experiments. 


Experiment 

Assimilation 

Cycle 

Initialization 

Type 

Remarks 

1  -  Control 

N/A 

♦AVMI 

No  assimilation— uses 

NMC  analyses 

2  -  12DI 

12  Hours 

*DVMI 

Assimilation  without 
asynoptic  GALE  data 

3  -  06DI 

6  Hours 

DVMI 

Assimilation  with  6-hourly 
GALE  data 

4  -  03DI 

3  Hours 

DVMI 

Assimilation  with  3-hourly 
GALE  data 

5-06AI 

6  Hours 

AVMI 

Assimilation  with  6-hourly 
GALE  data;  no  diabatic 

heating  in  model  initialization 


*AVMI  -  Adiabatic  vertical  mode  initialization 
*DVMI  -  Diabatic  vertical  mode  initialization 
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forecast  served  as  the  first-guess  field  for  the  initial  analysis  in  the  assimilation 
sequences.  In  experiment  2,  the  intermittent  data  assimilation  was  performed  from 
25/1 2Z  to  28/OOZ  with  a  12-h  update  cycle,  where  only  synoptic  data  entered  into  the 
analyses.  An  update  cycle  of  6  h  was  used  in  experiment  3.  GALE  data  valid  at  00, 06, 
12,  and  18Z  were  incorporated.  Like  experiment  2,  the  assimilation  is  performed  for  a  2 
1/2-day  period.  Due  to  extremely  limited  sounding  data  at  25/15Z  and  25/21Z,  the 
assimilation  run  in  experiment  4  began  at  26/OOZ,  using  the  26/OOZ  analysis  from 
experiment  3  as  its  initial  assimilated  state.  The  update  cycle  here  is  3  h  to  use  all  3- 
hourly  GALE  data  from  26/OOZ  to  28/OOZ.  Diabatic  vertical  mode  initialization  was  used 
in  experiments  2,  3,  and  4.  Experiment  5  was  identical  to  experiment  3  with  a  6-h  update 
cycle,  except  diabatic  heating  was  excluded  fi’om  the  vertical  mode  initialization.  The 
assimilated  states  (analyses  at  each  update  cycle)  from  experiments  2-5  were  saved  for 
post-processing,  i.e.  inter-comparisons  and  comparisons  with  the  NMC  analyses.  For 
all  experiments,  a  3-h  data  window  centered  on  the  analysis  time  was  used. 

From  the  assimilated  fields  valid  at  26/OOZ  and  26/1 2Z  of  each  experiment,  24-h 
forecasts  were  produced  using  the  same  model  that  was  used  in  the  assimilation.  Prior  to 
model  integration,  these  assimilated  analyses  were  initialized  using  the  vertical  mode 
initialization  procedure.  The  model  initialization  using  analyses  from  experiments  2, 3, 
and  4  included  forcing  due  to  diabatic  heating.  But  for  the  sake  of  fairness,  the  diabatic 
heating  was  derived  from  only  rainfall  occurring  prior  to  forecast  start  time  (a  3-h 
accumulation),  versus  a  6-h  accumulation  centered  on  the  analysis/initialization  time  as 
done  in  the  assimilation  runs. 
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6.4  Results 

Evaluation  of  the  experiments  is  performed  by  means  of  data-fitting  statistics, 
where  rms  differences  between  observations  and  analyses  and  between  observations  and 
background  forecast  fields  are  compared.  The  evaluation  also  includes  a  more  subjective 
comparison  of  synoptic- scale  and  mesoscale  patterns,  amplitudes,  and  positions  of 
specific  features  from  the  analyses  (assimilated  states)  and  the  subsequent  24-h  forecasts. 

6.4. 1  Statistical  Evaluation  and  Discussion 

Root-mean-square  errors  for  the  fit  of  background  forecasts  and  analyses  to  the 
observations  at  each  cycle  of  the  assimilation  experiments  were  compiled.  The  data  used 
here  are  those  passed  by  the  data  checking  procedures  in  the  analysis  scheme.  The  bogus 
data  are  not  included  in  this  statistical  evaluation.  These  statistics  were  computed  for  each 
analysis  variable  (geopotential,  u-  and  v-  wind  components,  temperature,  and  relative 
humidity)  at  each  of  the  19  analysis  levels  from  1000  mb  to  100  mb.  Figure  6.2  shows 
the  average  analysis  and  average  background  forecast  rms  errors  in  geopotential  for  the 
assimilation  analyses  and  forecasts  valid  at  0000  and  1200  UTC  from  experiments  2,  3, 
and  4  which  had  update  cycles  of  12, 6,  and  3  h,  respectively.  By  extracting  the  statistics 
at  only  the  synoptic  h''urs  (cycles),  the  comparison  is  more  equitable  than  if  including 
statistics  from  all  hours  of  experiments  3  and  4.  In  Fig.  6.2b,  the  rms  error  in  the 
background  geopotential  at  the  synoptic  hours  was  greatest  at  all  levels  when  using  the 
12-h  update  cycle.  In  the  upper  troposphere  and  at  the  surface,  the  3-h  update  cycle  led 
to  the  least  error  with  about  5  -  10  m  less  error  than  with  the  12-h  cycle.  These  results 
show  that  a  positive  impact  is  made  on  the  first  guess  field  of  geopotential  at  the  synoptic 
hours  when  3-hourly  GALE  data  are  included.  However,  in  the  mid-  and  lower 


Figure  6.2.  (a)  Average  analysis  rms  errors  and  (b)  average  background  forecast  rms 
errors  for  geopotendal  (m)  from  the  0000  and  1200  UTC  update  cycles  of  experiments  2 
(12DI),  3  (06DI).  and  4  (03DI). 
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troposphere,  the  6-h  update  cycle  gave  slightly  better  results  than  when  updating  every  3 
h.  A  possible  explanation  could  be  that  the  3-h  cycle  may  not  allow  sufficient  time  for 
model  spinup  in  regions  not  updated  by  the  3-hourly  asynoptic  data.  The  corresponding 
average  analysis  rms  errors  at  the  synoptic  hours,  displayed  in  Fig.  6.2a,  are  significantly 
less,  as  expected,  than  the  background  errors.  Here,  the  inclusion  of  the  3-hourly  data 
led  to  a  decrease  in  analysis  errors  throughout  the  troposphere,  with  the  assimilation 
using  the  3-h  update  cycle  being  better  than  the  6-h  update  cycle. 

Error  statistics  were  also  compiled  at  the  asynoptic  hours  of  0600  and  1800  UTC 
from  experiments  3  and  4.  As  shown  in  Fig.  6.3,  the  background  and  analysis  errors  in 
geopotential  decreased  at  all  levels  as  additional  3-hourly  data  were  assimilated;  that  is, 
the  assimilation  with  an  update  cycle  of  3  h  showed  an  improvement  over  the  case  with  a 
6-h  cycle.  Clearly  at  the  asynoptic  hours,  the  improved  background  fields  from  the  high- 
resolution  forecast  model  are  successful  in  carrying  information  from  earlier  asynoptic 
observations  forward  in  time,  leading  to  an  improved  analysis  which  then  results  in  a 
better  background  forecast  for  the  next  update  cycle  of  the  assimilation. 

Figures  6.4  -  6.6  illustrate  the  average  rms  statistics  of  u-,  v-,  and  relative 
humidity  from  experiments  2  -  4  for  the  synoptic  update  cycles.  Even  though  the 
analysis  error  in  the  wind  and  humidity  fields  remained  essentially  unchanged  as  the 
assimilation  cycle  decreased,  the  subsequent  forecasts  improved  at  all  levels  with  the 
decreasing  cycle  length. 

Average  root-mean-square  errors  for  geopotential  from  experiments  3  and  5  (all 
cycles  included)  are  plotted  in  Fig.  6.7.  Assimilation  with  diabatic  heating  included  in  the 
initialization  procedure  (06DI)  resulted  in  background  errors  within  the  lower  troposphere 
of  about  1  m  less  than  in  the  assimilation  excluding  the  convective  heating  (06AI),  but 
little  or  no  improvement  at  higher  levels.  The  analysis  errors  from  these  two  experiments 
were  almost  identical.  Similar  results  were  found  for  the  wind  components  and  relative 
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Figure  6.4,  Average  analysis  and  background  forecast  rms  errors  for  the  u-  wind 
component  (m  s"^)  from  the  0000  and  1200  UTC  update  cycles  of  experiments  2  (12DI), 
3  (0®I),  and  4  (03DI).  Analysis  error  plots  are  on  the  left  side  of  the  figure. 


component  (m  s‘l)  from  the  0000  and  1200  UTC  update  cycles  of  experiments  2  (12DT), 
3  (06DI),  and  4  (03DI).  Analysis  error  plots  arc  on  the  left  side  of  the  figure. 
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Figure  6.6.  Average  analysis  and  background  forecast  rms  errors  for  relative  humidity 
(%)  from  the  0000  and  1200  UTC  update  cycles  of  experiments  2  (12DI),  3  (06DI),  and 
4  (03DI).  Analysis  errw  plots  are  on  the  left  side  of  the  figure. 
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Figure  6.7.  (a)  Average  analysis  rms  errors  and  (b)  average  background  forecast  rms 
errors  for  geopotendal  (m)  from  all  update  cycles  of  experiments  3  (ObDI)  and  5  (06AI). 
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humidity.  These  results  suggest  that  the  inclusion  or  exclusion  of  the  diabatic  forcing  in 
the  assimilation  model  initialization  has  minimal  impact  on  the  assimilation  of  these  first 
order  variables. 

Generally,  the  background  forecast  errors  at  all  levels  (as  shown  in  Figs.  6.2  - 
6.7)  are  in  good  agreement  with  and  verify  our  estimates  of  the  standard  deviations  of  the 
observations  from  the  first  guess  forecasts  (refer  to  Table  6.3).  These  estimates  are  a 
sum  of  the  input  observational  error  and  a  linear  forecast  growth  error.  With  further 
testing,  we  can  utilize  statistics  compiled  over  a  longer  period  to  further  adjust  these  input 
error  statistics  for  use  in  the  analysis. 

Further  insight  can  be  gained  by  investigating  the  time  variation  of  the  statistics  in 
the  case  of  the  assimilation  with  update  cycle  of  3  h  (experiment  4).  Figure  6.8a  shows 
plots  of  average  rms  analysis  errors  for  the  geopotential  at  0300  and  1500  UTC,  0600 
and  1800  UTC,  0900  and  2100  UTC,  and  1200  and  0000  UTC  from  experiment  4.  The 
analysis  error  (differences  between  analysis  and  observations)  is  less  at  the  synoptic 
hours  of  0000  and  1200  UTC  than  at  the  asynoptic  hours.  The  larger  errors  noted  by 
Rodgers  et  al.  (1990)  in  the  geopotential  data  from  the  supplemental  GALE  soundings 
along  the  east  coast  may  have  contributed  to  the  higher  rms  errors  at  the  asynoptic  hours. 
At  the  synoptic  hours,  more  data  were  used  covering  the  whole  analysis  domain  (see  Fig. 
6.1),  and  the  average  analysis  error  would  be  less.  The  average  rms  background  errors 
for  gccjpotential  are  shown  in  Fig.  6.8b.  In  general,  the  background  fields  for  the  0000 
and  1200  UTC  analyses  had  the  largest  rms  error,  while  those  for  the  0300  and  1500 
UTC  analyses  possessed  the  least  error.  From  the  smaller  analysis  error  and  the  larger 
data  coverage  at  the  synoptic  hours,  one  would  expect  a  more  accurate  background 
forecast  for  the  0300  and  1500  UTC  analyses.  At  synoptic  hours,  the  background  in 
regions  outside  of  the  influence  of  the  asynoptic  data  would  be  less  accurate  and  the 
comparison  with  synoptic  data  would  lead  to  larger  errors.  However,  at  0600  and  18(X) 
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TABLE  6.3.  Standard  deviation  of  the  observational  errors  and  updated  forecast  error 
growth  for  a  6-h  forecast  period  in  parentheses.  These  are  the  input  error  statistics  for  the 
objective  analysis. 


Pressure 

Height 

Temperature 

Rel  humidity 

u,v 

Wind 

(mb) 

(m) 

(*C) 

(%) 

(m 

100 

20  (25) 

2.0  (0.7) 

3.5 

(1.3) 

150 

18 (25) 

1.8  (0.6) 

4.0 

(1.4) 

200 

15  (20) 

1.8  (0.5) 

4.0 

(1.5) 

250 

14(18) 

1.5  (0.6) 

4.0 

(1.6) 

300 

14(16) 

1.0  (0.6) 

17(4) 

4.0 

(1.6) 

350 

14(14) 

1.0  (0.5) 

15(4) 

4.0 

(1.7) 

400 

12  (12) 

1.0  (0.5) 

13(4) 

3.5 

(1.7) 

450 

10(11) 

1.0  (0.4) 

13(4) 

3.5 

(1.6) 

500 

9(10) 

1.0  (0.3) 

13(4) 

3.0 

(1.6) 

550 

8(10) 

1.0  (0.3) 

13(3) 

3.0 

(1.5) 

600 

7  (9) 

1.0  (0.3) 

12  (3) 

3.0 

(1.4) 

650 

6  (8) 

1.2  (0.3) 

12(3) 

3.0 

(1.3) 

700 

6  (8) 

1.3  (0.3) 

12(2) 

2.5 

(1.2) 

750 

6  (8) 

1.4  (0.3) 

12(2) 

2.5 

(1.2) 

800 

6  (8) 

1.5  (0.4) 

13(2) 

2.5 

(1.1) 

850 

6  (8) 

1.5  (0.5) 

13(2) 

2.5 

(1.1) 

900 

5  (8) 

1.8  (0.5) 

13(2) 

2.5 

(1.1) 

950 

5  (8) 

1.8  (0.6) 

13(2) 

2.5 

(1.1) 

1000 

5  (8) 

1.8  (0.6) 

13(2) 

2,5 

(1.1) 

Figure  6.8.  (a)  Average  analysis  rms  errors  and  (b)  average  background  forecast  rms 
errors  for  geopotential  (m)  from  experiment  4  (03DI)  at  the  0300  /  1500  UTC,  0600  / 
1800  UTC,  0900  /  2100  UTC,  and  1200  /  0000  UTC  update  cycles. 
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UTC,  the  errors  in  the  background  fields  were  larger  than  the  errors  in  the  first-guess 
fields  valid  at  0900  and  2100  UTC. 

Similar  statistics  were  compiled  from  experiment  4  for  the  u-  and  v-  wind 
components  and  relative  humidity  (not  shown).  However,  the  results  are  not  as  clear  cut 
as  in  the  case  of  geopotential.  The  background  error  at  the  0300  and  1500  UTC  cycles 
was  generally  less  than  that  at  the  synoptic  cycles,  and  likewise  the  analysis  error  at  the 
0300  and  1500  UTC  cycles  was  generally  larger  than  that  at  the  0000  and  1200  UTC 
update  cycles.  But  the  errors  of  the  two  sets  of  intermediate  cycles  (0600/1800  UTC  and 
0900/2100  UTC)  were  very  random,  with  values  not  always  falling  between  the  errors 
from  the  0300/1500  cycles  and  the  synoptic  cycles  as  in  the  case  with  geopotential. 

These  results  can  be  useful  in  adjusting  the  input  background  forecast  error  growth 
rate  used  in  the  objective  analysis  scheme.  Clurrently,  the  analysis  error  and  the  error 
growth  are  fixed  values  for  all  cycles  within  a  given  assimilation  sequence.  When 
considering  the  above  results,  it  may  be  prudent  in  the  future  to  adjust  the  analysis  error 
and  the  forecast  error  growth,  particularly  for  geopotential,  in  a  fashion  that  allows  the 
input  background  error  to  be  largest  at  the  0000  and  1200  UTC  update  cycles  and  least  at 
the  0300  and  1500  UTC  cycles. 

6.4.2  Subjective  Evaluation  and  Discussion 

Some  synoptic  analysis  of  the  results  is  beneficial  for  understanding  the  impact  of 
four-dimensional  data  assimilation  on  synoptic  and  mesoscale  features,  i.e.  in  the  present 
study,  the  impact  of  assimilating  the  3-hourly  GALE  observations.  Here,  emphasis  is 
placed  on  a  few  significant  events  in  the  surface  and  500-mb  analyses  and  the  forecast 
precipitation  fields.  Precipitation  forecasts  are  verified  against  our  observed  precipitation 
(raingauge  /  satellite-derived  rainfall).  Precipitation  verification  is  especially  valuable 
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since  rainfall  is  the  result  of  many  complex  physical  processes  and  is  usually 
characterized  by  considerable  small-scale  variability.  Since  rainfall  is  dependent  on 
divergence  and  vertical  motion,  it  is  also  a  good  indicator  of  the  accuracy  of  the  model's 
dynamic  balance  and  intervariable  consistency. 

a)  Analysis  aspects 

Sea-level  pressure  and  lOOO-mb  wind  and  temperature  analysis  fields  from  the  "no 
assimilation"  case  and  the  assimilation  with  a  12-h  update  cycle  (12DI)  at  25/12Z  are 
shown  in  Fig.  6.9.  The  12DI  analysis,  shown  in  Fig.  6.9b,  captured  the  strong 
temperature  gradient  associated  with  the  coastal  front  along  the  Carolina  coast  and  also 
the  cold  air  damming  with  strong  ridging  east  of  the  Appalachians.  The  25/1 2Z  Bosart 
surface  analysis  (Fig.  6.10)  verifies  the  existence  of  these  mesoscale  features.  However, 
the  "no  assimilation"  analysis,  shown  in  Fig.  6.9a,  has  a  weaker  thermal  gradient  along 
the  coast.  Strong  onshore  flow,  which  is  readily  apparent  in  Bosart's  analysis,  was  also 
present  in  the  assimilated  analysis.  Figure  6.11  shows  the  1000-mb  specific  humidity 
fields  from  experiments  1  and  2.  The  low-level  moisture  gradients  at  the  coastal  front 
and  with  the  short  wave  over  the  western  Gulf  of  Mexico  are  much  greater  in  the  12DI 
analysis  than  in  the  "no  assimilation"  analysis. 

At  26/OOZ,  surface  and  500-mb  analyses  from  experiments  1,  2,  3,  and  5  were 
examined.  Compared  to  the  "no  assimilation"  analyses,  the  assimilated  surface  fields 
possessed  much  stronger  thermal  and  low-level  moisture  gradients  across  the  coastal 
front,  with  little  distinction  between  the  three  assimilations  of  varying  update  cycles  (not 
shown).  The  500-mb  relative  humidity  analyses,  shown  in  Fig.  6.12,  reveal  a  better 
distinction.  While  rain  was  falling  over  central  North  Carolina  at  26A)0Z,  an  upper-level 
moisture  maxima  was  present  there  in  the  assimilated  fields.  However,  the  analysis  from 


Figure  6.10.  GALE  andysis  of  sca-lcvel  pressure  deviation  (ftom  1000  mb),  surface 
tempei^m  (  C)  and  wind  for  25/12Z  by  Bosart  (1988).  Surface  temperature  contours 
every  D 
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the  6-h  update  cycle  had  the  greatest  value  showing  that  asynoptic  GALE  data,  in  this 
case,  led  to  an  improvement  over  the  assimilation  with  12-h  updates.  Over  the  central 
Gulf  coast,  all  three  assimilations  produced  a  moisture  tongue  (with  relative  humidity 
over  70%)  associated  with  an  eastward-moving  short-wave  disturbance.  The  "no 
assimilation"  analysis  showed  a  maxima  of  only  56%  in  this  region  where  rainfall  was 
occurring.  The  overall  upper-level  moisture  pattern  from  the  06AI  assimilation  (i.e., 
VMI  without  diabatic  heating)  was  very  close  to  the  other  assimilated  states.  As  noted  in 
the  previous  section,  the  inclusion  of  diabatic  heating  in  the  vertical  mode  initialization  led 
to  little  change  in  the  first  order  field.  This  is  not  necessarily  true  in  the  case  of  vertical 
motion,  which  is  related  directly  to  diabatic  heating.  From  the  vertical  velocity  fields, 
shown  in  Fig.  6.13,  a  clear  distinction  between  the  experiment  with  diabatic  VMI 
(experiment  3)  and  that  with  adiabatic  VMI  (experiment  5)  can  be  seen.  The  assimilation 
with  diabatic  heating  included  in  VMI  shows  a  stronger  vertical  circulation  at  the  coastal 
front  at  26/OOZ.  The  vertical  motion  along  the  east  coast  from  the  assimilated 
experiments  was  much  stronger  and  had  much  more  detail  than  that  from  the  "no 
assimilation"  case,  with  the  06DI  assimilation  producing  slightly  greater  vertical  motions 
at  the  coastal  front  than  did  the  assimilation  with  a  12-h  update  cycle.  Also,  stronger 
vertical  velocities  are  present  over  the  central  Gulf  coast  associated  with  the  developing 
wave  in  all  the  assimilated  states  as  compared  to  the  "no  assimilation"  analysis. 

Cross  sections  were  produced  to  further  investigate  the  ageostrophic  circulation 
associated  with  the  coastal  front  and  an  advancing  upper-level  jet.  Cross  section  2  (CS2 
from  90* W,  43.5*N  to  68.7° W,  26°N)  cuts  across  the  entrance  region  of  the  jet  over  the 
southeast  U.S.  and  extends  across  the  Carolina  coast  (see  Fig.  4.5c).  Figure  6.14  shows 
the  magnitude  of  the  wind  speed  normal  to  the  cross  section  from  the  "no  assimilation" 
case  and  the  12DI  assimilation  for  26/OOZ.  The  ageostrophic  circulations  in  the  plane  of 
CS2  for  these  two  cases  are  depicted  in  Fig.  6.15.  The  mesoscale  circulation  ot  the 


Figure  6.13.  500-mb  vertic^  velocity  (mb  h*l)  analyses  at  26fQ07.  %)m  (a)  experiment  1 
(no  assimilation),  (b)  experiment  2  (12DI),  (c)  experiment  3  (06DI),  and  (d)  experiment  5 
(06AI).  Vertical  motion  of  magnitude  of  mb  h'^  or  more  are  contoui«i  every  2.5  mb 
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Figure  6.14.  Wind  speed  (m  s'l)  normal  to  the  plane  of  CS2  at  26/OOZ  from  (a) 
expenment  1  (no  assimilation)  and  (b)  experiment  2  (12DI).  Contours  are  every  5  m  s*^. 
The  potential  temperature  (*K)  is  also  shown  by  the  dashed  contours  every  5  *K. 
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Hgure  6.15.  Vertical  motion  (mb  h*l)  and  ageostrc^hic  wind  (m  s'l)  in  the  plane  of  CS2 
at  26/OOZ  from  (a)  experiment  1  (no  assimilation)  and  (b)  experiment  2  (12DI).  The 
magnitude  of  the  vertical  motion  is  represented  by  the  solid  contours  at  -40,  -25,  -20,-15, 
-12.5  and  every  2  mb  h'l  for  -10  and  above. 
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coastal  front  is  still  a  distinct  feature  within  the  12DI  analysis.  However,  the  jet’s 
secondary  circulation  has  begun  to  absorb  and  deepen  the  coastal  front's  ageostrophic 
circulation  as  shown  in  the  cross  section.  The  ageostrophic  circulation  at  26/OOZ  due  to 
the  coastal  front  is  much  stronger  and  sharper  in  the  12DI  as  compared  to  the  control 
analysis,  with  maximum  ascents  of  -18.2  and  -8.7  mb  h‘^  at  the  coastal  front, 
respectively.  The  06DI  assimilation  (cross  section  not  shown)  strengthened  the 
ageostrophic  circulation,  but  only  by  a  small  amount  (-19.1  mb  h'^  vs.  -18.2  mb  h'l). 

At  26/1 2Z,  there  is  little  difference  in  the  surface  temperature,  wind,  and  pressure 
fields  from  the  various  assimilations  sequences.  In  the  1000-mb  specific  humidity 
fields,  the  moisture  gradient  along  the  front  on  the  east  coast  was  sharpest  in  the  3-h  cycle 
assimilation,  followed  by  06DI,  then  12DI,  with  the  "no  assimilation"  analysis  having  the 
weakest  gradient.  The  500-mb  relative  humidity  analyses  from  experiments  1-4  again 
depicted  the  tongue  of  upper-level  moisture  extending  over  the  eastern  U.S.  from 
Georgia  to  Maine.  As  in  the  lower  troposphere,  the  03DI  assimilation  generated  the 
sharpest  gradient  and  largest  maxima  and  the  "no  assimilation"  case  produced  the 
weakest.  Little  difference  was  found  in  the  vertical  motion  fields  resulting  from  the 
assimilation  runs,  with  all  displaying  much  greater  maxima  and  detail  than  in  the  "no 
assimilation"  experiment. 

At  27/C)0Z  and  27/1 2Z,  the  assimilated  fields  from  experiments  2, 3,  and  4  are  very 
similar.  The  influence  of  the  3-hourly  GALE  data  seems  to  be  impacting  the  assimilation 
primarily  when  mesoscale  circulations  are  predominant,  i.e.,  there  is  a  positive  impact  in 
the  small-scale.  When  the  synoptic-scale  circulation  became  stronger  and  essentially 
engulfed  the  coastal  front,  less  impact  was  noted.  However,  each  of  the  assimilations 
still  produced  analyses  which  were  noticeably  different  from  the  "no  assimilation"  case. 
One  surface  feature  worth  noting  at  27/OOZ  is  the  positions  of  the  low  pressure  centers 
from  each  experiment.  By  21/OOZ  the  disturbance,  which  had  moved  from  the  northern 


Gulf  of  Mexico  across  southern  Georgia,  again  appeared  as  a  closed  surface  low  in  the 
assimilation  analyses;  rapid  cyclogenesis  had  just  begun  (Fig.  6.16a).  The  low  pressure 
center  which  was  earlier  associated  with  the  coastal  front  had  moved  north  and  would 
soon  merge  with  a  weak  northern  low,  associated  with  a  short-wave  trough  which  had 
advanced  eastward  from  the  Great  Lakes  region  into  New  England.  However,  the  "no 
assimilation"  analysis  does  not  distinguish  those  circulations,  instead  only  the  deepening 
low  pressure  center  is  depicted  in  a  position  farther  up  the  east  coast  (Fig.  6.16b).  By 
27/03Z,  the  03DI  assimilation  had  deepened  the  cyclone's  central  pressure  by  3  mb  (Fig. 
6.16c)  and  moved  it  just  off  the  Delaware  coast  in  close  proximity  with  the  "no 
assimilation"  position  at  27/OOZ.  This  would  seem  to  suggest  that  only  a  slight  timing 
error  in  the  development  of  the  strengthening  cyclone  differentiates  the  "no  assimilation" 
analysis  from  the  assimilation  case.  However,  investigation  of  the  ageostrophic 
circulations  along  the  east  coast  proves  otherwise.  As  in  other  cases  of  rapid 
cyclogenesis  (Uccellini  and  Kocin  1987,  for  example),  the  rapid  development  is  being 
aided  by  the  indirect  circulations  of  an  exiting  upper-level  jet  to  the  north  and  an 
approaching  jet  located  south  of  the  low.  In  this  case,  the  northern  jet  has  also  been 
strengthened  by  the  short-wave  trough  from  the  (jreat  Lakes.  These  jets  are  evident  in 
Fig.  6.17a  which  is  a  cross  section  (CS7,  see  Fig.  6.16a)  along  the  east  coast  v69.3*W, 
49.5'’N  to  80.0’W,  26.5*N)  of  the  normal  wind  speed  at  27/OOZ.  Figures  6.17b  and 
6.17c  are  cross  sections  of  the  ageostrophic  circulations  from  the  03DI  and  "no 
assimilation"  cases.  The  03DI  analysis  shows  the  ascending  motion  associated  with  the 
nonhem  low  (point  A  of  Fig.  6.17b),  still  retains  the  mesoscale  circulation  associated 
with  the  decaying  coastal  low  (point  B),  and  depicts  the  strong  vertical  circulation  of  the 
developing  low  at  point  C.  The  secondary  circulation  of  the  coastal  low  is  being 
absorbed  by  the  indirect  circulation  of  the  entrance  region  of  the  northern  jet.  The 
ageostrophic  circulation  from  the  "no  assimilation"  analysis  (Fig.  6.17c)  fails  to  capture 
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Figure  6.16.  Sea-level  pressure  (mb),  1000-mb  temperature  (*Q  and  wind  analyses  at 
27/OOZ  from  (a)  experiment  4  (03DI)  and  (b)  experiment  1  (no  assimilation),  and  at 
27A)3Z  from  (c)  experiment  4  (03DI).  Contours  as  in  Fig.  6.9. 
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Figure  6.17.  (a)  Cross  section  (CS7)  of  the  normal  wind  speed  (m  s'^)  at  27/OOZ  from 
experiment  4  (03DI).  CS7  of  the  ageostiophic  circulation  at  27/OOZ  firom  (b)  experiment 
4  (03DI)  and  (c)  experiment  1  (no  assimilation)  and  at  27/03Z  from  (d)  experiment  4. 
Contours  as  in  Figs.  6.14  and  6.15.  See  Fig.  6.16a  for  location  of  CS7. 
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these  various  features.  It  simply  depicts  the  strong  vertical  motion  of  the  deepening  low, 
being  supported  by  the  indirect  circulations  of  the  two  jets.  Figure  6.17d  shows  the  03DI 
ageostrophic  circulation  at  27/03Z.  By  this  time,  the  analyzed  deepening  low  had  moved 
a  little  farther  north  (point  C).  The  circulations  associated  with  the  northern  low  and 
coastal  low  are  merging  and  continue  to  be  engulfed  by  the  jet's  secondary  circulation. 

By  27/1 2Z,  the  rapidly  deepening  cyclone  had  tracked  northeast  to  New  England 
in  both  the  "no  assimilation"  and  03DI  cases  (Figs.  6.18a  and  6.18b).  With  the  suppon 
of  the  secondary  circulations  of  the  two  jets,  the  low  had  deepened  by  15  mb  in  the  "no 
assimilation"  analysis  and  17  mb  in  the  03DI  analysis.  The  synoptic  scale  circulation  of 
the  cyclone  was  now  totally  dominant,  with  the  mesoscale  circulation  of  the  coastal  low 
being  completely  absorbed  in  the  large-scale  circulation.  As  a  result,  the  pattern  of  the 
ageostrophic  circulation  along  CS6  (68.7‘W,  55.5*N  to  76.7*W,  25.5*N;  see  Fig. 
4.24b)  from  the  03DI  assimilation  (Fig.  6.19b)  is  very  similar  to  that  from  the  "no 
assimilation"  case  (Fig.  6.19a).  But  the  cyclone's  vertical  motion  from  the  assimilation 
case  is  much  stronger,  with  a  maximum  of  -45  mb  h‘l  versus  -26.2  mb  h’^. 

b)  Subsequent  forecast  aspects 

The  true  test  of  the  goodness  of  any  data  assimilation  designed  for  the  purpose  of 
providing  initial  conditions  for  a  numerical  weather  prediction  model  is  the  accuracy  of 
the  forecasts  produced  from  the  assimilated  states.  Twenty-four-hour  forecasts  starting  at 
26/OOZ  were  generated  from  analyses  in  experiments  1  (no  assimilation),  2  (12DI),  3 
(06DI),  and  5  (06AI).  Comparisons  among  these  forecasts  were  made  and  are  now 
discussed. 

Six-hour  forecasts,  valid  at  26/06Z,  of  sea-level  pressure,  surface  air  temperature 
and  winds  are  shown  in  Fig.  6.20.  The  forecasts  originating  from  the  assimilated  fields 
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Figure  6.18.  Sea-level  pressure  (mb),  1000-mb  temperature  (*C)  and  wind  analyses  at 
IIUTL  from  (a)  experiment  1  (no  assimilation)  and  (b)  experiment  4  (03DI).  Contours 
as  in  Fig.  6.9. 
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Figure  6.19.  Cross  sections  (CS6)  of  the  ageostrophic  circulation  at  27/1 2Z  from  (a) 
experiment  1  (no  assimilation)  and  (b)  experiment  4  (03DI)-  Contours  as  in  Fig.  6.15. 
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Figure  6.20.  Six-hour  forecasts  of  sea-level  pressure  (mb),  surface  air  tempmture  (*Q, 
and  1000-mb  winds  valid  at  26A)6Z  from  (a)  experiment  1  (no  assimilation),  (b) 
experiment  2  (12DI),  and  (c)  experiment  3  (06DI).  Contours  and  in  Fig.  6.9. 
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have  much  tighter  temperature  gradients  along  the  Carolina  coast  in  the  vicinity  of  the 
coastal  front  than  does  the  forecast  from  the  "no  assimilation"  case.  The  frontal 
temperature  gradients  resulting  from  experiments  2  and  3  are  in  very  good  agreement 
with  the  Bosart  hand  analysis  (Fig.  6.21)  with  a  temperature  contrast  of  approximately 
15“C  from  Charlotte,  North  Carolina  to  a  position  just  off  the  South  Carolina  coast. 
Bosart's  analysis  shows  a  low  of  1012  mb  on  the  Florida  panhandle  and  a  low  of  less 
than  1012  mb  which  had  formed  on  the  southern  end  of  the  coastal  front.  In  reference  to 
these  low  pressure  centers,  the  12DI  forecast  verifies  the  best  in  regards  to  central 
pressure,  but  the  06DI  forecast  predicts  the  positions  more  accurately.  As  in  the  case  of 
the  analyses,  the  forecasts  from  the  assimilated  states  produce  a  much  tighter  low  level 
moisture  gradient  at  the  coastal  front  than  does  the  "no  assimilation"  forecast  (not 
shown). 

Six-hour  observed  rainfall  and  accumulated  forecast  rainfall  (26/OOZ  -  26/06Z)  for 
experiments  1,3,  and  5  are  shown  in  Fig.  6.22.  The  forecasts  from  the  assimilated 
cases,  06DI  (Fig.  6.22b)  and  12DI  (not  shown),  match  the  observed  precipitation  pattern 
(Fig.  6.22d)  reasonably  well  and  show  much  more  rapid  spinup  of  rainfall  than  does  the 
"no  assimilation"  forecast.  Here,  the  disadvantage  of  not  having  a  mesoscale  assimilation 
model  to  generate  and  maintain  small-scale  features  is  evident.  The  "no  assimilation" 
forecast  must  spinup  these  circulations  and,  as  a  result,  gives  a  poor  short-range  rainfall 
prediction.  The  forecast  from  experiment  5  (06AI  in  Fig.  6.22c)  fails  to  capture  any 
rainfall  over  the  Gulf  coast  associated  with  the  developing  short  wave.  The  absence  of 
diabatic  heating  from  the  VMI  has  caused  a  spinup  problem  in  this  case. 

The  12-h  surface  and  low-level  moisture  forecasts  (valid  at  26/1 2Z)  from  the 
various  assimilations  are  very  similar  and,  like  the  6-h  forecasts,  have  much  stronger 
gradients  along  the  east  coast  of  the  U.S.  than  does  the  "no  assimilation"  forecast. 
Figure  6.23  depicts  the  6-h  accumulated  rainfall  ending  at  26/1 2Z.  The  forecasts  from 


Figure  6.21.  GALE  andysis  of  sea'level  pressure  deviation  (from  ICKX)  mb),  surface 
temperature  (*C)  and  wind  for  26/06Z  by  Bosart  (1988).  Surface  temperature  contours 
every  5’C. 
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Figure  6.22.  Forecasts  of  accumulated  total  precipitation  in  cm  (26AX)Z  -  26A)6Z)  from 
(a)  experiment  1  (no  assimilatitm),  (b)  experiment  3  (06DI),  and  (c)  expoiment  5  (06AI). 
(d)  Observed  precipitation  in  cm  (26/OOZ  -  26A)6Z).  Contours  at  0.25  cm  and  every  0.5 
cm  for  0.5  and  ab^e. 
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06DI  and  12DI  are  much  better  than  that  from  the  "no  assimilation"  case  and  verify 
reasonably  well  with  the  observed  rainfall.  However,  the  forecast  from  12DI  produced  a 
maximum  of  only  1.5  cm  near  Cape  Hatteras,  whereas  the  06DI  forecast  maximum  was 
nearly  3  cm,  a  value  much  closer  to  the  observed  maximum  of  5.31  cm.  The  rainfall 
forecast  from  experiment  5  (06AI)  continued  to  be  too  weak  on  the  southern  extent  of  the 
precipitation  area  and  the  "no  assimilation"  case  almost  totally  missed  this  feature.  The 
forecasts  from  all  the  experiments  completely  missed  the  maximum  of  6.85  cm  in  New 
Jersey  as  depicted  in  Fig.  6.23d.  The  horizontal  resolution  of  even  our  fine-resolution 
mesoscale  model  is  insufficient  to  resolve  individual  convective  cells.  The  model's 
cumulus  parameterization  can  not  be  expected  to  duplicate  the  observed  precipitation 
maxima  resulting  from  these  sub-grid  scale  convective  storms. 

By  26/1 8Z  the  accumulated  6-h  precipitation  forecasts  (Fig.  6.24)  had  become 
more  similar  with  each  extending  the  rainfall  area  across  central  Florida.  However,  the 
()6DI  and  12DI  forecasts  did  extent  the  precipitation  area  over  most  of  North  Carolina  and 
South  Carolina  in  agreement  with  the  observed  pattern,  whereas  the  "no  assimilation" 
forecast  had  moved  all  precipitation  offshore  for  latitudes  south  of  Virginia. 

A  second  set  of  forecasts  with  an  initiation  time  of  26/1 2Z  were  produced  from 
experiments  1  -  4.  Figure  6.25  depicts  6-h  forecasts  of  sea-level  pressure  and  lOOO-mb 
temperature  and  wind  fields;  these  forecasts,  including  the  one  from  the  "no  assimilation" 
analysis,  are  very  similar.  Note  similarities  in  temperature  gradients  along  the  east  coast 
of  U.S.  and  the  locations  of  the  coastal  low  defrned  by  the  center  of  the  low  pressure 
over  southeast  New  York  state.  The  coastal  low  is  seen  merging  with  a  trough  of  low 
pressure  to  the  northwest.  The  new  cyclone  center  is  forming  near  Cape  Hatteras  in  the 
extended  trough  along  the  coast.  The  low-level  moisture  fields  (not  shown)  are  equally 
similar,  panicularly  the  sharp  gradient  associated  with  the  frontal  system  along  the  east 
coast.  The  6-h  accumulated  rainfall  forecasts  ending  at  26/1 8Z  originating  from  the 


Hgure  6.24.  Forecasts  of  accumulated  total  precipitation  in  cm  {26I12L  -  26/18Z)  from 
(a)  experiment  1  (no  assimilation),  (b)  expoimnt  2  (12DI),  and  (c)  experiment  3  (06DI). 
(d)  Observed  precipitation  in  cm  (26/12Z  >  26/1 8Z).  Contours  same  as  in  Fig.  6.22. 
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Figure  6.25.  Six-hour  forecasts  of  sea-level  pressure  (mb)  and  10(X)-mb  temperature 
(*C)  and  wind  valid  at  26/1 8Z  firom  (a)  experiment  1  (no  assimilation),  (b)  experiment  2 
(12DI),  (c)  experiment  3  (06DI),  and  (d)  experiment  4  (03DI). 
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assimilated  fields  are  much  better  than  the  forecast  originating  firom  the  "no  assimilation" 
analysis,  with  the  forecasts  from  the  3-h  (03DI,  not  shown)  and  6-h  (06DI)  cycle 
assimilations  being  the  closest  to  and  in  good  agreement  with  the  observed  rainfall  pattern 
(Fig.  6.26).  In  other  words,  the  assimilation  process  had  an  impact:  there  is  a  definite 
improvement  over  the  "no  assimilation"  case,  and  the  assimilation  of  asynoptic  3-hourly 
GALE  data  had  at  least  a  minor  impact  in  the  rainfall  since  the  assimilations  with  the  3-h 
and  6-h  cycles  led  to  slightly  better  forecast  precipitation  fields  than  the  case  from  the  12- 
h  cycle  assimilation. 

During  the  remainder  of  the  prediction  period,  the  synoptic-scale  circulation  of  this 
developing  east  coast  cyclone  was  dominant  and  the  forecast  patterns  from  the  different 
experiments  were  quite  similar.  The  12-  and  24-h  forecast  precipitation  fields  from  each 
case,  including  the  "no  assimilation"  experiment,  were  similar  and  in  reasonable 
agreement  with  the  observed  precipitation.  By  27/OOZ,  a  low  pressure  center  appeared 
off  the  North  Carolina  coast  and  deepened  rapidly  as  it  moved  along  the  coast  (Fig. 
6.27).  In  the  12-h  forecasts  (valid  at  27/()0Z)  from  the  assimilated  fields,  the  region  of 
low  pressure  extends  farther  north  and  indicates  two  different  northern  centers  (Fig. 
6.27b),  which  compare  more  favorably  with  the  assimilated  analysis  shown  in  Fig. 
6.16a.  The  24-h  surface  forecasts  from  the  control  and  03DI  cases  are  shown  in  Figs. 
6.27c  and  d.  These  forecasts  are  very  similar,  except  that  the  center  of  the  low  has  a 
better  position  in  the  assimilated  case  when  compared  with  the  analysis  in  Fig.  6.18b.  It 
seems  that  the  assimilation  of  the  GALE  data  had  less  impact  particularly  in  the 
precipitation  during  the  latter  part  of  the  forecast  period.  The  rainfall  was  being  forced 
primarily  by  the  large  synoptic-scale  circulation  of  the  developing  east  coast  cyclone, 
enabling  the  control  forecast  to  produce  adequate  results.  However,  the  12-  and  24-h 
surface  pressure  forecasts,  with  the  assimilation  of  the  GALE  data,  show  more 


Figure  6.26.  Forecasts  of  accumulated  total  precipitation  in  cm  (26/1 2Z  -  26/1 8Z)  from 
(a)  experiment  1  (no  assimilation),  (b)  expeiiment  2  (12DI),  and  (c)  experiment  3  (06DI), 
with  forecast  start  time  of  26/12^  (d)  Observed  precipitation  in  cm  (26/12Z  -  26/1 8Z). 
Contours  same  as  in  Fig.  6.22. 


Figure  6.27.  12-h  forecasts  of  sea-level  pressure  (mb)  and  lOOO-mb  temperature  (*C) 
and  wind  valid  at  27A)0Z  ftom  (a)  experimnt  1  (no  assimilation)  and  (b)  experiment  4 
(03DI),  and  24-h  forecasts  valid  at  27/12Z  from  (c)  experiment  1  (no  assimilation)  and 
(d)  experiment  4  (03DI). 
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consistency  with  the  assimilated  analyses  than  the  control  forecasts  do  with  the  NMC 
analyses. 
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7.  CONCLUSIONS 

Four-dimensional  data-assimilation  methodology  has  been  applied  to  high- 
frequency  observational  data  using  a  multivariate,  successive  correction  objective  analysis 
scheme,  a  diabatic  vertical  mode  initialization  procedure,  and  a  fine-resolution  mesoscale 
numerical  weather  prediction  model.  With  only  a  limited  number  of  iterations,  the 
NRL/NCSU  objective  analysis  scheme,  based  on  Bratseth  weighting  functions, 
converges  to  the  optimal  interpolation  solution.  This  scheme  produces  accurate  analysis 
fields  by  using  wind  information  to  achieve  gradients  of  geopotential  through  the 
geostrophic  relationship.  It  is  also  successful  in  correcting  the  phase  of  the  first-guess 
temperature  field,  improving  the  second-order  fields  of  absolute  vorticity  and  vertical 
motion,  strengthening  the  first-guess  upper-level  jet,  and  resolving  mesoscale 
circulations. 

When  diabatic  vertical  mode  initialization  is  performed  on  these  analysis  fields, 
model  spinup  is  reduced.  With  the  inclusion  of  diabatic  heating,  which  forces  convective 
precipitation,  the  VMI  retains  the  divergence,  and  therefore  vertical  motion,  associated 
with  the  convection.  This  diabatic  initialization  procedure  is  used  to  initialize  the 
background  forecasts  at  each  cycle  of  the  assimilation  process.  Observed  rainfall, 
coupled  with  "spunup"  model  rainfall,  is  used  to  derive  three-dimensional  heating  fields. 
The  method  ensures  temporal  consistency  of  the  diabatic  heating  rates  by  using  rainfall 
rates  which  are  time-centered  (-3H  to  -»-3H)  on  the  analysis/initialization  valid  time.  This 
procedure  also  ensures  that  the  heating  rates  in  the  initial  integration  of  the  subsequent 
assimilated  background  forecast  are  representative  of  the  "observed"  heating. 

The  performance  of  the  assimilation  system  was  evaluated  by  means  of  data-fitting 
statistics  and  subjective  comparisons  of  assimilated  states  and  subsequent  forecasts.  The 
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assimilated  background  forecast  rms  error  for  all  analysis  variables  decreased 
significantly  as  the  length  of  the  update  cycle  shortened,  i.e.  as  additional  asynoptic, 
high-density  GALE  data  were  assimilated.  The  same  was  true  for  the  analysis  error  of 
geopotential.  This  decrease  in  background  and  analysis  error  clearly  shows  that  the 
inclusion  of  the  3-hourly  GALE  data  had  a  positive  impact  on  the  assimilation.  The 
analysis  errors  of  the  wind  components  and  relative  humidity  varied  little  between  the 
different  experiments,  due  most  likely  to  relatively  large  observational  error  in  the  wind 
and  humidity.  The  background  and  analysis  rms  errors  changed  little  when  diabatic 
heating  was  included  in  the  vertical  mode  initialization,  suggesting  that  the  inclusion  or 
exclusion  has  minimal  impact  on  the  assimilation  of  these  first-order  variables.  The 
assimilation  with  a  3-h  update  cycle  revealed  that  the  background  forecast  error  generally 
increases  as  the  period  from  the  previous  synoptic  update  to  the  ciurent  analysis/forecast 
cycle  increases.  With  only  small  insertions  of  data  at  the  asynoptic  cycles  (as  compared 
to  complete  datasets  at  the  synoptic  hours),  the  effectual  background  forecast  duration  of 
successive  update  cycles  increases  between  synoptic  updates  for  geopotential. 

A  subjective  evaluation  of  the  assimilation  system  was  conducted  by  comparing 
surface  and  500-mb  analyses  and  forecast  fields.  The  analyses  and  subsequent  forecasts 
from  each  of  the  assimilation  experiments  were  consistently  better  than  the  "no 
assimilation"  analyses  (NMC  analyses)  and  subsequent  forecasts.  A  temporal 
consistency  was  achieved  between  the  assimilated  analyses  and  forecasts,  which  was  not 
achieved  without  assimilation.  Stronger  thermal  and  moisture  gradients  and  vertical 
motions  were  evident  at  the  coastal  front  in  the  assimilated  fields.  There  was  less 
difference  between  the  various  assimilations  which  included  the  diabatic  VMI.  The 
surface  analyses  of  pressure,  wind,  and  temperature  and  vertical  motion  fields  were 
generally  very  similar.  However,  the  C3DI  and  06DI  lower-  and  upper-level  moisture 
analyses  showed  some  improvement  over  the  12DI  analyses.  Here,  the  inclusion  of  the 


206 


3-hourly  GALE  data  made  a  positive  impact  on  the  moisture  analyses.  Excluding  diabatic 
heating  from  the  VMI  (06AI)  led  to  no  noticeable  detriment  in  the  analyses  of  first  order 
variables,  but  did  result  in  weaker  vertical  velocities  associated  with  the  mesoscale 
circulation  of  the  coastal  front  and  upper-level  jet  streaks.  This  result  is  reasonable  since 
vertical  motion  is  directly  related  to  diabatic  heating. 

A  subjective  verification  of  forecast  6-h  accumulated  precipitation  was  conducted. 
Precipitation  verification  is  a  valuable  measure  of  the  forecast  "goodness"  since  rainfall  is 
dependent  on  many  complex  physical  processes;  as  a  result,  it  is  much  more  difficult  to 
predict  accurately  than  the  other  modelled  variables.  The  forecast  originating  from  the 
assimilated  states  verified  reasonably  well  with  the  observed  rainfall  and  were  much 
better  than  the  "no  assimilation"  forecasts,  showing  much  faster  spinup  of  forecast 
precipitation  than  did  the  control  case.  The  assimilations  with  shorter  update  cycles 
(included  asynoptic  GALE  data)  led  to  slightly  improved  precipitation  forecast  as 
compared  to  the  12DI  forecasts. 

The  results  of  this  study  show  that  the  assimilation  of  GALE  data  in  a  mesoscale 
model  is  successful  in  producing  significantly  improved  analyses  and  subsequent 
forecasts  which  are  more  consistent  in  time.  The  results  also  indicate  that  a  small  area  of 
higher  temporal  data  resolution  can,  when  assimilated  with  high-resolution  mesoscale 
prediction  model,  lead  to  improved  assimilated  states  and  better  forecasts.  However,  the 
impact  of  the  GALE  data  on  the  forecasts  was  generally  greater  when  mesoscale 
circulations  were  important,  such  as  in  the  interaction  of  the  coastal  front  and  the  upper- 
level  jet  during  the  period  prior  to  rapid  cyclogenesis.  When  the  large-scale  synoptic 
system  became  dominant,  the  3-hourly  data  had  less  impact  in  the  subsequent  forecasts. 

In  the  present  study,  only  diabatic  heating  associated  with  convective  precipitation 
is  considered.  Future  work  should  include  the  investigation  of  latent  heating  resulting 
from  large-scale,  non-convective  rainfall.  Also,  forecast  precipitation  totals  in  areas  of 
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strong  convection  were  sometimes  far  from  the  observed  amounts  even  when  the  3- 
hourly  GALE  data  were  assimilated.  The  extremely  small  scale  of  the  convective  activity 
precluded  our  model  from  resolving  these  feamres.  Work  has  begun  to  increase  the 
model's  horizontal  and  vertical  resolutions  and  to  include  vertical  accelerations.  That  is, 
the  NRL/NCSU  model  will  be  converted  to  a  16-layer,  non-hydrostatic  model  with 
horizontal  grid  spacing  on  the  order  of  5  to  10  km.  Significant  improvement  in  rainfall 
forecasts  should  result. 
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